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laminates  displayed  conplex  distortions  of  their  initial  yield  surface  during  plastic  work-hardening.  Analysis  of  a 
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SECTION  1 


INTRODUCTION 


When  studying  the  deformation  of  metallic  solids  due  to  applied 
external  loads,  it  is  common  to  assume  the  material  will  behave  in  an 
isotropic  manner  [  1  ].  Typically,  this  assumption  is  extended  to  both 
the  elastic  and  e1astic>p1astic  response  [2].  In  many  cases,  the 
particular  manufacturing  process,  such  as  rolling,  introduces 
anisotropy  into  the  material,  giving  directionally  sensitive  elasto- 
plastic  mechanical  behavior  [3]. 

There  are  two  categories  of  problems  involving  anisotropic 
behavior.  One  is  the  induction  of  anisotropy  into  an  initially 
isotropic  bulk  material  [4  ],  most  coninonly  observed  during  metal 
forming  operations.  The  second  category  involves  those  materials 
which  are  initially  anisotropic.  An  example  is  the  behavior  of  some 
metal -matrix  composites.  Both  categories  of  problems  are  pertinent 
to  this  research  and  a  finite  element  program  ANPLAST  is  developed  to 
this  end.  This  program  implements  a  properly  formulated  anisotropic 
plasticity  theory  in  a  convenient  form  to  use  for  numerical  purposes. 

The  anisotropic  elastoplasticity  of  sheet  materials  is  often 
studied  from  two  perspectives.  The  first  takes  into  consideration  the 
transverse  or  through- thickness  anisotropy  while  assuming  planar 
•  isotropy.  Most  of  the  research  to  date  has  concentrated  on  this 

subject  [3,5]  because  this  is  the  usual  case  incurred  in  rolling. 

While  the  problem  of  planar  anisotropy  has  been  investigated  from  a 
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theoretical  perspective  [4,6,7],  very  few  examples  of  its  applications 
can  be  found  in  the  literature.  Consequently,  this  research  will 
concentrate  exclusively  on  planar  anisotropy  and  its  influence  on 
some  practical  problems  of  interest. 

Specifically,  a  thin  sheet  with  a  circular  hole  and  a  panel  with 
a  crack  will  be  studied.  The  differences  between  the  isotropic  and 
anisotropic  solutions  are  elaborated.  The  elastoplastic  analysis  is 
extended  to  layered  (or  laminated)  structures  with  a  particular 
application  to  the  sheet  with  a  circular  hole.  The  present  numerical 
solutions  of  these  problems  are  confirmed  by  comparison  with  experimen¬ 
tal  results.  Some  comparisons  are  also  made  with  the  available 
theoretical  solutions. 
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SECTION  2  FUNDAMENTALS  OF  ORTHOTROPIC  PLASTICITY 


INTRODUCTION 

This  section  will  review  the  fundamentals  of  elasticity  of 
orthotropic  materials.  A  similar  approach  will  then  be  used  to  derive 
the  orthotropic  plasticity  flow  rule.  The  meaning  of  the  coefficients 
that  appear  in  the  theory  will  be  demonstrated  for  uniaxial  response 
of  the  material.  A  discussion  of  the  path  dependency  associated  with 
plastic  flow  necessarily  entails  a  comparison  of  incremental 
deformation  theories  of  plasticity  and  this  is  also  covered  in  this 
section.  In  anticipation  of  the  later  applications  to  laminated 
composites,  formulation  of  the  elastoplastic  response  of  a  multilayered 
laminate  is  included. 

BASICS  OF  ELASTICITY 

For  an  elastic  body,  the  current  stress  state  depends  only  on  the 
current  deformation  state.  This  can  be  expressed  as  [  1  ] 


aw 


(2.1) 


where 

=  stress  tensor 
“  strain  tensor 
W  *  energy  density  function 
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For  linear  elasticity.  W  is  assumed  quadratic  in  the  stress  as 


“  *  ^ijki  °ij 
Equation  2.1  then  becomes 

®ij  “  ^ijki  °ki 


(2.2) 


The  tensor  of  elastic  constants  is  independent  of  stress  and 

strain  and  exhibits  the  following  symmetric  behavior 


^ijki  *  ^jikt'  ^ijkt  *  ^ijik'  ^ijki  *  ^ktij 


The  matrix 


has  a  maximum  of  21  independent  constants. 


which  can  be  written  in  contracted  notation  as 


If  the  material  has  three  mutually  perpendicular  planes  of  syametry, 
it  is  considered  orthotropic  and  the  stress-strain  relation  reduces  to 


4 


having  9  independent  constants  as. 


When  analyzing  thin  sheets  of  material,  the  through-thickness 
stress  can  be  neglected  as  well  as  the  corresponding  shear  stresses, 
resulting  in  a  state  of  plane  stress.  The  constitutive  equation  now 
becomes 


A  total  of  4  independent  constants  describe  the  planar  problem. 
The  through  thickness  strain  e^^  could  be  obtained  from  knowledge  of 
and  C22.  The  coefficients  in  Eq.  2.3  can  be  represented  in  the 
more  common  engineering  form  as 

^11  “  -q  ^12 ' '  if 
hz  •  q  Se  ”  ^ 
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where 


E^,  E2  *  elastic  moduli  1 
V12  *  Poisson  ratio 

G^2  ~  shear  modulus 

and  subscripts  1,  2  denote  the  corresponding  principal  material  directions. 

Since  the  material  coefficients  are  constants,  the  incremental 
constitutive  relation  for  the  plane  orthotropic  case  simply  becomes 


where  [c]  has  the  same  meaning  as  Eg.  2*3. 

The  incremental  load  deformation  response  of  an  elastic  material 
is  completely  defined  by  Eg.  2.4. 

PLASTICITY  FORMULATION 

The  approach  taken  to  the  plasticity  formulation  is  to  use  a  flow 
rule  that  is  guadratic  in  the  stresses.  This  gives  sufficient 
flexibility  to  describe  orthotropic  plasticity  L  2  ]. 

Three  Dimensional  Anisotropic  Flow  Rule 

The  incremental  plastic  strains  can  be  written  in  terms  of  a 
plastic  potential  [ 8  ]  as 


6 


where  the  superscript  "p"  denotes  plasticity.  For  an  anisotropic  body, 
Kachanov  [8  ]  suggested  that  g  take  the  quadratic  form 


It  is  interesting  to  note  that  a  number  of  authors  [9,10,11  ]  have 
investigated  non-quadratic  potential  functions,  but  a  quadratic 
function  is  deemed  sufficient  for  the  purposes  of  this  research. 
Taking  advantage  of  the  symmetric  properties  of  Eq.  2.5  is 

contracted  and  rewritten  (as  in  the  elastic  case)  as 


'“'ll 

"Cn 

Ci2  Cj3  ....  Cjg 

r  's 
°n 

“'22 

C22 

°Z2 

"'33 

>  = 

i 

^33! 

"'23 

Symmetric 

\ 

°23 

"'13 

°13 

"'12 

^66 

°12 

The  matrix  [C^j]^  contains  21  independent  coefficients.  If  the 
material  contains  three  mutually  perpendicular  planes  of  symmetry, 
then  the  plastic  coefficient  matrix  takes  the  form 


[C]” 


^11 

^12 

^13 

0 

0 

0 

^12 

C22 

^23 

0 

0 

0 

^13 

^23 

^33 

0 

0 

0 

0 

0 

0 

^44 

0 

0 

0 

0 

0 

0 

^55 

0 

0 

0 

0 

0 

0 

^66 

Thus,  to  describe  three  dimensional  orthotropic  plastic  flow  requires 
only  nine  independent  coefficients. 


Reduction  to  Two  Dimensional  Plane  Stress 

When  dealing  with  plane  stress,  the  plastic  flow  relations 
simplify  further  since  0^3  *  0,  0^3  *  0,  033  ®  0,  giving 


“Cu  C12  o~ 

r 

/  de22 

>  ■ 

C12  C22  0 

”22 

i?®12^ 

0  0  ^66 

! 

°12 

(2.7) 


The  above  system  of  equations  are  over-specified  with  four  C-values  and 
dA*.  Normalize  these  equations  by  introducing  the  following 
definitions 

r  =  /C^ 

*"66  *  ^66^^22  *  ^22 

r^Z  •  ^12^^22 
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Then,  the  incremental  expressions  for  plane  stress  become 


‘"ll'j 

a. 

■"ii 

''12 

o 

°ll'| 

/de22 

• 

^12 

1 

0 

< 

>°22  ^ 

de,, 

0 

0 

''66 

b2> 

(2.8) 


It  is  clear  from  this  equation  that  the  incremental  plastic  strains  are 
orthotropic  with  respect  to  the  stresses,  i.e.,  the  normal  components 
of  incremental  strain  are  uncoupled  from  the  shear  stress.  This 
relation  reduces  to  isotropic  plastic  flow  when  the  r's  have  the 
values  [  12  ] 


rii  =  1.0  rj2  =  -.50  rgg  *  1.50 


The  factor  of  proportionality  dA  is  solved  from  work  considerations. 
Parenthetically,  the  plane  stress  orthotropic  flow  rule  could  have  been 
derived  directly  by  taking  the  potential  function  g  as 


g  =  I  [r 


11®11 


+  a 


22 


^'"66°12 


^'"l2°ll‘^22^ 


Definition  of  Effective  Stress  and  Strain 


In  general,  the  plastic  work  increment  can  be  written  as 


which  for  plane  stress  expands  to 


dW^  =  ®22^®22  *  ^°12*^®12 


The  plastic  strain  increments  previously  defined  in  Eq.  2.8  can  be 
substituted  into  the  above  relation  resulting  in 


dW^  =  *  °22  *  ^*^66^12  *  ^'^12°11°22^  *  2gdA 


Let  the  effective  stress  and  strain  be  defined  such  that 


diP  =  dwP  X  2gdA 


Further,  the  effective  stress  can  be  defined  as 


2  -2 
3  " 


^'■66°12  *  ^’‘l2°ll®22  * 


Therefore , 


dx  *  ^  (diP/a)  *  ^  (deP/do)(do/a) 


and  the  effective  plastic  strain  can  be  written  as  [12] 


[diP]  »  I  [defi  >  r^jdePg  -  2rj2dePideP2 


2(rii  -  r^2)  ”  '^12  ^ 


(2.9) 


(2.10) 


(2.11) 


(2.12) 


(2.13) 
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Substituting  the  above  expression  for  dx  into  Eq.  2.8,  the  plane  stress 
orthotropic  plastic  flow  equations  become 


""11 

'"12 

1 - 

o 

d®22 

• 

*"12 

1 

0 

0 

0 

"'66. 

do  0 

(2.14) 


No  assumption  of  constant  volume  deformation  or  incompressibility 
was  made  during  the  development  of  these  plastic  flow  equations.  To 
consider  volumetric  changes  or  thickness  changes  in  the  2-D  case 
would  require  additional  coefficients. 

Layered  Elastoplastic  Materials 

A  laminate  is  defined  as  two  or  more  laminas  bonded  together  to 
act  as  an  integral  structural  element.  Let  the  laminate  consist  of 
thin  layers,  each  exhibiting  plane  stress  behavior.  Due  to  the 
potential  for  plastic  flow  in  each  layer,  all  stress-strain  relations 
will  be  written  in  incremental  form.  The  elastic  and  plastic  incremen¬ 
tal  constitutive  relations  can  be  written  in  matrix  form  as 

=  [Q]®  {de®},^ 

(2.15) 

{do}|^  *  [Q]^ 

where  [Q]®  and  [Q]^  are  the  transformed  inverse  of  matrices  of  Eqns. 

2.4  and  2.14,  respectively  for  the  layer  of  a  multilayered  laminate. 
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The  total  stress-strain  relation  Is  obtained  by  adding  the  elastic 
and  plastic  responses  to  give 

{do}^  =  [Q],^  {de,^}  (2.16) 

where 


[Q]  =  CQ]'  ♦  [Q]" 

It  Should  be  noted  that  [Q]  will  be  stress  dependent  when  plastic 
deformation  occurs  through  the  contribution  of  [Q]^. 

Knowledge  of  the  variation  of  stress  and  strain  through  the 
laminate  thickness  is  essential  to  the  extensional  stiffness  of  a 
laminate.  For  the  purposes  of  this  research,  bending  stiffness  and 
deformation  will  be  Ignored  and  only  the  first  order  extensional 
behavior  considered.  The  laminate  Is  assumed  to  consist  of  perfectly 
bonded  lamlnas  with  displacements  continuous  across  the  lamina 
boundaries  so  that  no  lamina  can  slip  relative  to  another.  These  are 
the  usual  type  of  assumptions  made  In  laminate  theory  [13]. 

A  typical  laminate  Is  pictured  in  Fig.  2.1  along  with  Its  deformed 
shape.  It  Is  quite  clear  that  the  deformation  and  strain  In  the  plane 
of  each  layer  is  the  same.  The  incremental  stress-strain  relations 
for  each  lamina  can  be  expressed  as 


{da}^  .  [Q]^  {de^,} 


(2.17) 
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Figure  2.1  Illustration  of  an  N-layered  laminate  along 
with  a  typical  deformed  geometry. 
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where  {de^}  is  the  laminate  strain  increment. 

The  resultant  laminate  forces  per  unit  width  are  obtained  by 
integrating  the  stress  components  of  each  lamina  through  the  total 
thickness  of  the  laminate.  They  will  take  the  following  form 


/  > 

^  \ 

do 

XX 

(  dN 

)  = 

i 

do„„ 

y 

] 

-t/2 

yy 

dN 

do„.. 

L 

Since  the  stresses  are  constant  through  the  thickness  of  each  lamina, 
then  this  becomes 
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Replacing  the  incremental  lamina  stresses  with  those  from  Eq.  2.17, 
results  in  the  incremental  laminate  load-deformation  equations  taking 
the  form 
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(2.18) 
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The  global  response  of  the  laminate  will  be  affected  by  each 
layers  stiffness  contribution  [Q]|^-  During  elastoplastic  deformation, 
each  layer  accounts  for  its  particular  plastic  stiffness  through  [Q]^- 

ORTHOTROPIC  BEHAVIOR  UNDER  UNIAXIAL  STRESS 
The  effect  of  the  plasticity  coefficients  appearing  in  the 
previous  constitutive  relations  can  best  be  understood  by  analyzing  a 
uniaxial  stress  state. 

Interpretation  of  the  Plasticity  Coefficients 

For  a  uniaxial  specimen  cut  from  an  orthotropic  material  at  an 
arbitrary  direction  to  the  principle  axes  as  illustrated  in  Fig.  2.2,  the 
state  of  stress  with  respect  to  the  material  directions  are  as 
follows 


®11 

°22 

°12 


OgCos^e 

.  2 
OgSin  e 

OgSinecose 


The  incremental  strains  in  the  loading  direction  and  in  the  transverse 
direction  are 

de^  *  [r^jcos^e  +  2rj2COS^8sin^e  +  sin^e  +  2rggSin^ecos^e]  dx 
de^  «  [r^j  +  l-2rgg)  sin^ecos^e  +  r^g  (sin^e  +  cos^e]  dx 

(2.19) 
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Figure  2.2  Illustration  of  a  uniaxial  specimen  cut  at 
an  arbitrary  angle  to  the  principal  material 
directions. 


The  ratio  of  transverse  to  longitudinal  strain  increments  can  be 
viewed  as  a  plastic  Poisson  ratio  and  is  written 


p  *  l-2rgg)sin^ecos^e  +  rj2(sin^e  +  cos^e) 

®  dej  r^jcos^e  +  sin^e  +  2(rj2  +  rgg)sin^ecos^e 

(2.20) 

Note  that  for  various  orientations  the  sine  and  cosine  terms  are  zero, 
and  the  above  expression  simplifies,  i.e. 
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e  *  90° 


e  «  45°  :  -vP  , 
45‘ 


'"12 


(2.21) 


Thus,  it  seems  that  the  plasticity  coefficients  behave  similarly  to 
Poisson  ratio's.  If  the  isotropic  values  of  the  r's  are  substituted 
into  the  above  equation,  the  plastic  Poisson  ratio  is  ^  for  all 
orientations.  Otherwise,  if  the  plasticity  coefficients  are  varied 
from  their  isotropic  values,  anisotropy  will  be  induced  into  the 
plastic  flow. 

Parametric  Study  of  the  Plasticity  Coefficients 

Additional  insight  can  be  obtained  by  conducting  a  simple 
parametric  study  of  the  three  plasticity  coefficients.  The  axial  and 
transverse  loading-deformation  behavior  of  an  arbitrarily  oriented 
uniaxial  specimen  will  be  studied.  The  elastic  behavior  of  the 
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material  will  be  assumed  to  be  isotropic  so  that  the  anisotropy  of  the 
plastic  flow  can  be  highlighted. 

A  FORTRAN  computer  program  MULTIAX  was  developed  for  the  purpose 
of  predicting  the  load-deformation  behavior  of  an  ideal  element  under 
multiaxial  loading.  This  program  calculates  the  incremental  strains 
during  each  load  increment  and  sums  them  to  obtain  the  total  strain. 

The  program  can  handle  any  loading  path. 

The  plasticity  parameter  r^j  is  studied  first.  Figure  2.3 
illustrates  the  load-deformation  behavior  with  assumed  isotropic 
properties.  As  expected,  the  global  stress-strain  response  is  not  a 
function  of  e.  The  isotropic  example  will  be  used  as  the  base  response 
for  comparison.  Changing  the  value  of  r^^  to  .8  and  then  .4  as  shown 
in  Fig.  2.4,  highlights  the  global  stress-strain  dependence  on  e.  It 
is  quite  evident  that  the  more  r^^^  deviates  from  unity,  the  stronger 
the  induced  anisotropy.  It  is  obvious  that  yield  does  not  occur  at  a 
common  value  for  all  orientations.  This  is  because  the  yield  is  a 
function  of  the  r's. 

Figure  2.5  depicts  the  uniaxial  load-deformation  as  r^2  changed 
from  its  isotropic  value  of  -0.5.  Figure  2.6  shows  the  effect  as  rgg 
is  changed  from  its  isotropic  value  of  2.5.  The  departure  of  either 
of  these  plasticity  coefficients  from  their  isotropic  value  also 
induces  anisotropy  into  the  uniaxial  yield. 

These  studies  show  that  each  parameter  can  affect  the  anisotropy 
of  material  behavior. 

Interelation  of  Work-hardening  and  the  Plasticity  Coefficients 

Plasticity  coefficients  play  an  important  role  in  determining 
the  unique  relationship  between  effective  stress  (o)  and  the  effective 
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Figure  2.3  Uniaxial  elastoplastlc  response 
of  an  isotropic  tnaterial. 
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Figure  2.4  Uniaxial  elastoplastic  response  of  a  material 
in  which  the  plasticity  coefficient  r,.  is 
varied. 


Figure  2.4  continued 
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Figure  2,5  Uniaxial  elastoplastic  response  of  a  materia! 
in  which  the  plasticity  coefficient  r,,  is 
varied. 
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Figure  2.5  continued 


i 
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Figure  2.6  Uniaxial  elastoplastic  response  of  a  material 
in  which  the  plasticity  coefficient  rgg  is 
varied. 
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plastic  strain  (i^).  These  coefficients  allow  the  stress-strain  curves 
obtained  at  different  angles  to  be  collapsed  into  a  single  master  curve. 
With  this  in  mind,  the  work-hardening  relationship  for  any  arbitrary 
orientation  of  a  uniaxial  specimen  can  be  written  as 


0  =  h(i^) 


(2.22) 


where 


°  *  's[\  ^‘’a 

'VF  i  ‘e 

V  4  4  2  2  2  2 

r^jcos  0  +  sin  e  +  2rggSin  ecos  6  +  2rj2S'’»i  ®cos  6 


For  a  material  with  a  given  anisotropy  (e.g.  particular  r's),  the 
workhardening  is  specified  by  the  function  h(e^).  This  function  is 
usually  expressed  in  the  form  of  a  power  law.  An  example  of  such  a 
hardening  law  is  a  three  parameter  model  used  by  Kenaga  [12]  to 
represent  a  Boron-Aluminum  metal -matrix  composite.  It  appears  as 


eP 


(2.23) 


where  o^,  a,  and  6  are  curve  fitting  parameters. 

In  the  work  to  follow,  the  r's  are  constants;  so  the  function 
h(iP)  is  the  only  means  by  which  the  work  hardening  can  be  changed. 
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It  is  interesting  to  note  in  Figs.  2. 3-2.6  that  all  parametric  cases 
utilize  an  identical  work -hardening  function,  yet  the  individual 
uniaxial  elastic-plastic  slopes  are  different. 


YIELD  CRITERION  UNDER  MULTIAXIAL  LOADING 
For  the  case  of  an  ideally  isotropic  material  under  a  uniaxial 
stress  state,  a  single  value  of  yield  will  suffice.  But  a  multiaxial 
stress  complicates  the  definition  of  the  yield  point,  and  it  actually 
becomes  a  yield  surface.  A  further  complication  is  the  introduction  of 
anisotropy  into  the  material.  It  is  apparent  that  the  locus  of  yield 
points  under  various  combinations  of  stress  forms  a  surface  in  stress 
space. 

A  possible  interpretation  of  effective  stress  (o)  is  a  single 
surface  in  space  for  a  single  point  on  the  work-hardening  curve. 
Therefore,  one  can  say  that  yield  will  occur  at  a  value  of  effective 
stress,  say  o^,  and  is  a  function  of  the  plastic  strain.  That  is,  the 
general  yield  criterion  can  be  written  as 
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Oil  ^  On'i  +  2r,oOi,  Ooo  +  2r,,o"  ) 
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66'’12 


(2.24) 


In  general,  this  expression  defines  the  yield  surface  for  an 
orthotropic  material  under  a  multiaxial  state  of  stress. 

Parametric  Study  of  the  Plasticity  Coefficients 

In  order  to  provide  a  graphical  representation  of  the  orthotropic 
ellipsoid,  the  yield  surface  is  plotted  in  022*  0^,2  space.  A 
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parametric  study  of  r^^,  illustrated  in  Fig.  2.7. 

The  study  of  r^j  indicates  a  definite  distortion  of  the  ellipsoid 
in  the  direction  of  increasing  o^j  for  decreasing  r^^.  It  is  clear 
that  r^2  (which  multiplies  the  stress  cross  term  0^1022^  affects  the 
symmetric  elongation  along  the  plane  =  022*  l^his  result  implies 
that  yield  in  a  uniaxial  stress  state  oriented  along  a  principle 
material  axis  will  not  be  affected  by  rj2*  But  yield  in  a  multiaxial 
stress  will  be . 

The  plasticity  coefficient  rgg  can  be  studied  by  observing  the 
influence  of  0^2*  ^  decreasing  r^g  elongates  the 

ellipsoid  in  the  direction  of  the  shear  stress  component.  This  has 
the  effect  of  decreasing  the  contribution  of  0^2  towards  yielding. 

Isotropic  Behavior 

It  is  of  interest  to  consider  the  requirements  necessary  for 
isotropic  yield  behavior.  A  material  is  isotropic  if  the  yield 
criterion  is  invariant  under  a  rotation  of  axes.  That  is 


-2 


■  const 


3  2^2 

2  (>"11011  •>•  022 

3  '2  *  ‘2 
1  ®22 


+  2rj2°ll°22  ^’"66°12^ 
+  2rj2®li®22  ^’’66^12^ 


where  ojj,  022*  0^2  a>"®  In  the  rotated  system.  In  order  for  this  to 
be  true  for  any  transformed  stresses,  then  it  is  necessary  that 


rjj  .  1.0 


•■iZ  *  ''66  ■  ‘-O 


(2.25) 
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Figure  2.7  Illustration  of  the  yield  surface  in 

°11’  '^22*  ®12  space  for  different 

values  of  the  plasticity  coefficients. 
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It  is  interesting  to  note  that  this  is  more  general  than  the  usual 
isotropic  restrictions  of 

rii  *  1.0  r^2  “  -'50  '‘66  " 

This  point  is  demonstrated  in  Fig.  2.8  where  the  values 

'*11  *  *"12  * 

are  used.  These  curves  are  identical  to  those  of  Fig.  2,3. 

Effect  of  Work-hardening  on  the  Yield  Surface 
The  anisotropic  yield  locus  was  shown  to  be 

5^(5")  •  I  (r„cfi  ♦  0I2  *  '■l2’’U“22  * 

After  the  onset  of  initial  yielding,  plastic  deformation  will  occur. 

If  isotropic  hardening  is  assumed,  the  yield  surface  in  stress  space 
will  expand  to  each  new  stress  as  work>hardening  continues  as 
illustrated  in  Fig.  2.9.  Isotropic  hardening  results  in  equal  tension 
and  compression  for  a  uniaxial  specimen.  Note  that  the  yield  criterion 
is  anisotropic  even  though  the  hardening  is  isotropic.  This  is  often 
a  common  point  of  misunderstanding. 

Anisotropic  work-hardening  theory  accounts  for  changes  in  the 
shape  of  the  yield  locus  during  plastic  deformation.  Anisotropic 
hardening  can  be  introduced  by  letting  the  plasticity  coefficients 
vary  as  a  function  of  some  state  variable.  This  can  be  expressed  as 
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Figure  2.9  Illustration  of  isotropic  hardening  of 
an  anisotropic  yield  surface. 
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Oy(e^)  *  F(o, 

where  i  is  some  function  of  the  plasticity,  say  i^.  A  graphic 
depiction  of  this  hardening  behavior  is  presented  in  Fig.  2.10,  where 
not  only  is  changing,  but  so  are  r^j  and  r^2‘ 

LOAD  PATH  DEPENDENCE  OF  ACCUMULATED  STRAIN 
The  plasticity  theory  formulated  is  incremental.  The  advantage 
of  this  approach  is  that  it  allows  modelling  of  complex  loading  paths 
(including  unloading).  Hence,  a  discussion  of  path  dependence  is 
pertinent. 

Deformation  Theory 

The  deformation  theory  of  plasticity  developed  in  1924  by  Hencky 
assumes  a  unique  stress-strain  relationship.  As  presented  by  Kachanov 
[8] 


®ij 


ij 


^eP. 
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(2.26) 


where 


c  *  unit  volume  change  due  to  mean  pressure 
Tj  ■  effective  yield  stress 

and  r  Is  related  to  work-hardening.  For  example,  r  can  be  written  as 
a  function  of  the  effective  stress  or  effective  plastic  strain  as 
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0^=25000 
*11  =  -8 
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Figure  2^10  Illustration  of  anisotropic  hardening  of 
an  initially  isotropic  yield  surface. 
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r  *  f(o,  e^) 


It  was  also  shown  by  Kachanov  that  the  governing  equations  of 
deformation  theory  correspond  to  a  non-linear  elastic  constitutive 
representation.  Having  the  load-deformation  path  uniquely  defined 
results  in  significant  advantages  during  numerical  analyses,  in  that 
a  Newton-Rhapson  solution  algorithm  is  easily  implemented.  Unfortu¬ 
nately,  there  are  limitations  to  the  use  of  a  deformation  theory. 

Kachanov  demonstrated  that  the  governing  equations  of  deformation 
theory  coincide  with  incremental  theory  only  during  simple  loading. 

The  two  theories  agree  only  when  the  stress  components  maintain  the 
same  proportion  to  each  other.  A  simple  unloading  path  such  as  uniaxial 
unloading  would  violate  the  requirements  for  deformation  theory. 

Unloading  and  Cyclic  Loading 

A  monotonic  uniaxial  loading  and  corresponding  unloading  are 
depicted  in  Fig.  2.11.  One  can  imagine  that  the  loading  processes 
can  be  viewed  as  a  non-linear  elastic  material  response.  But  upon 
unloading,  it  becomes  quite  evident  that  the  elastoplastic  response 
is  not  well  represented  by  a  single  valued  stress-strain  relationship. 
Since  the  study  of  problems  with  elastic  unloading  and  residual 
stresses  is  of  primary  interest  in  this  research,  then  only 
incremental  theory  will  be  used. 

The  uniaxial  stress-strain  behavior  of  a  cyclic  loading  path  is 
pictured  in  Fig.  2.12.  This  illustration  further  reinforces  the  idea 
of  a  non-unique  stress-strain  response  and  the  limitations  of  the 
deformation  theory.  Parenthetically,  this  figure  also  highlights 
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Figure  2.11  Uniaxial  elastoplastic  response  of  an  isotropic 
material  which  is  loaded,  unloaded,  then 
reloaded. 
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Figure 
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.12  Uniaxial  elastoplastic  response  of  an  isotropic 
material  under  cyclic  loading. 
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the  effects  of  isotropic  hardening,  that  is,  the  self-similar  expansion 
of  the  yield  surface. 

This  parametric  study  highlights  the  complexity  of  the  yield 
criterion  in  a  multiaxial  stress  state.  If  anisotropy  were  introduced 
along  with  a  complex  geometry,  one  could  see  the  difficulty  in 
predicting  yield,  as  well  as  subsequent  work-hardening.  Because  of 
this  behavior.  Section  3 will  deal  with  the  implementation  of  the 
theory  into  a  finite  element  program. 

Some  Particular  Complex  Loading  Paths 

Complex  loading  paths  can  be  induced  by  the  application  of 
multiaxial  loads  in  sequence;  for  example,  a  uniaxial  load  followed 
by  a  transverse  load.  Figure  2.13  compares  three  different  load  paths 
involving  the  normal  stress  components.  Figure  2 . 14  on  the  other  hand, 
illustrates  the  effect  of  shear  along  with  a  normal  stress  component. 

In  both  figures,  the  case  of  simultaneous  load  application  (which 
would  be  the  deformation  case)  is  included.  These  results  clearly 
indicate  the  strong  path  dependence  of  accumulated  strain.  The  final 
stress-strain  response  is  obviously  not  unique  when  different  loading  paths 
are  taken.  This  rules  out  deformation  theory  for  this  class  of 
problems. 
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Figure  2.13  Elastoplastic  deformation  of  an  isotropic 
material  for  three  complex  loading  paths 
involving  both  normal  stress  components. 


Figure  2.14  Elastoplastic  deformation  response  of  an  isotropic 
material  for  three  complex  loading  paths  involving 
a  normal  and  a  shear  stress  component. 
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SECTION  3  FORMULATION  OF  COMPUTATIONAL  TOOLS 

INTRODUCTION 

The  study  and  application  of  orthotropic  plasticity  under  multi - 
axial  loading  requires  the  development  of  computational  tools.  Without 
these  tools,  it  is  very  difficult  and  expensive  to  gain  insight  into 
these  intricate  non-linear  problems.  The  incremental  formulation  of 
plasticity  makes  the  digital  computer  an  ideal  vehicle  to  investigate 
load-deformation  behavior  under  complex  loading  conditions. 

This  section  discusses  the  development  of  a  finite  element  program 
designed  to  analyze  orthotropic  elastic-plastic  problems.  The 
associated  pre-processing  and  post-processing  software  are  also 
developed.  The  lamination  theory  discussed  in  Section  2  is  extended 
and  included  in  a  computer  program  which  predicts  load  deformation 
behavior  of  laminated  orthotropic  elastic-plastic  materials. 

OVERALL  VIEW  OF  PROGRAM  ANPLAST 

A  self-contained  finite  element  program  and  peripheral  software 
written  in  Fortran  77  is  presented. 

Program  Structure  and  Philosophy 

The  ANPLAST  finite  element  program  is  designed  for  plane  stress 
anisotropic  elastic-plastic  analysis.  It  solves  the  non-linear  problem 
by  dividing  the  total  applied  load  into  small  load  increments  and  then 
summing  up  all  pertinent  quantities  such  as  stress,  strain. 
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displacement,  and  nodal  loads  [14].  The  global  stiffness  matrix  is 
updated  at  the  beginning  of  every  loading  increment  to  account  for  the 
inherent  material  non-linearity  associated  with  plasticity.  The 
analysis  continues  until  the  maximum  load  specified  is  reached,  or  as 
an  alternative,  the  number  of  yielded  elements  has  exceeded  a  maximum 
specified  by  the  user.  The  element  and  nodal  quantities  can  be  output 
in  a  number  of  optional  formats. 

A  descriptive  FORTRAN  flow  diagram  of  the  ANPLAST  program  is  shown 
in  Fig.  3  .1.  This  program  is  designed  so  that  the  main  program 
orovides  the  control  and  accounting  of  the  analysis  while  the  subroutines 
carry  out  the  main  functions.  The  list  of  subroutines  and  a  brief 
description  of  each  follows: 

INPUT:  Read  in  model,  geometry,  boundary  conditions 

and  initialize  the  program  control  parameters. 

MESG:  Provide  warning  and  error  messages. 

AMAT:  Form  the  element  strain-nodal  displacement  matrices. 

LOAD:  Read  in  the  incremental  nodal  loads. 

FORMK:  Form  the  global  stiffness  matrix. 

ESMAT:  Form  the  element  stiffness  matrix. 

SSE:  Form  the  incremental  elastic  stress-strain  relations. 

SSP:  Form  the  incremental  elastic-plastic  stress-strain 

relations. 

R0T4:  Angular  transformation  of  4th  order  tensor. 

R0T2:  Angular  transformation  of  2nd  order  tensor. 

INVERT:  Inversion  of  3x3  matrix 

OCMPBD:  Decomposition  of  global  stiffness  [15]  matrix. 

SLVBO:  Solution  of  global  system  of  equations  [15]. 
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Figure  3.1  Flow  Diagram  of  the  ANPLAST  finite  element  program. 
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STRESS 


Calculate  incremental  element  stresses  and  strains. 

OUTELM:  Output  element  stresses  and  strains. 

OUTNOO:  Output  nodal  stresses  and  strains. 

RESTRT:  Read  and  write  element/nodal  quantities  for  restart. 

Comparison  of  Solution  Algorithms 

Using  the  finite  element  technique  to  analyze  solids  results  in 
the  formation  of  a  large  number  of  simultaneous  equations.  There  is  a 
choice  to  include  the  solution  algorithm  directly  into  the  FORTRAN 
code,  or  utilize  scientific  library  utility  routines  (such  as  IMSL). 

It  was  decided  to  incorporate  the  equation  solver  directly  into  the 
ANPLAST  structure.  The  ANPLAST  program  consequently  is  self  contained 
and  therefore  does  not  require  any  peripheral  software  other  than  the 
usual  FORTRAN  77  libraries. 

A  parallel  version  of  ANPLAST  which  utilizes  the  IMSL  routines 
was  also  developed  to  verify  the  efficiency  of  the  solution  algorithm 
taken  from  Segerlind  [15].  The  computation  time  on  identical  analyses 
demonstrated  that  the  IMSL  version  was  actually  less  efficient  than  the 
hardwired  solution  routines.  Of  course  to  be  fair,  it  must  be 
realized  that  the  IMSL  routines  are  nrnre  general  in  nature,  and  easy  to 
access.  In  summary,  the  ANPLAST  solution  algorithms  are  efficient  and 
accurate  and  are  more  than  adequate  for  the  purposes  of  this  research. 

Program  Re»start  Capabilities 

The  majority  of  analyses  conducted  with  ANPLAST  involve  an  incre¬ 
mental  approach  to  the  non-linear  problem.  This  approach  usually 
requires  at  least  10  to  15  increments  and  results  in  considerable 
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computation  time.  In  addition,  the  size  of  each  loading  increment  must 
be  estimated  prior  to  the  analysis  and  may  turn  out  to  be  inappropriate. 
Further,  the  amount  of  data  which  can  be  printed  out  for  each  increment 
is  significant  and  usually  not  known  prior  to  the  analysis. 

Because  of  these  difficulties,  a  "restart*'  option  was  built  into 
ANPLAST.  This  "restart"  capability  permits  the  analysis  to  be  stopped 
and  started  again  at  a  later  time.  The  analysis  can  be  restarted  at 
whatever  increment  has  been  previously  written  to  the  restart  file. 

Two  blocks  of  information  are  stored:  the  element  quantities  and  the 
nodal  values.  The  analysis  is  re-started  by  reading  these  quantities 
in  and  re-initializing  the  appropriate  arrays.  New  incremental  loading 
is  supplied  and  the  analysis  continued. 

An  additional  benefit  of  the  restart  capability  is  the  post 
processing  of  the  output.  Only  a  minimal  amount  of  information  needs 
to  be  printed  out  after  each  increment,  since  the  analysis  can  be 
restarted  at  any  previous  increment  and  all  of  the  desired  information 
recovered. 

Numerical  Form  of  the  Incremental  Constitutive  Relations 

The  incremental  plastic  strains  represented  in  Eq.  2.14  have  to 
be  written  in  a  more  explicit  form  so  that  they  are  compatible  with  the 
finite  element  method.  Starting  with  the  equation  for  the  effective 
stress 


4 


I,  2 


+  0 


22 


^'"12®11°22 


^'"66®  12^ 


gives  the  effective  stress  increment  as 
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(3.1) 


do  =  w  -  [('"iiOn  +  r.^o-^)  do,,  +  (r. 


+  o,«)  do. 


2  "  "  ^12°22^  “°11  ''12”11  “22'  ““22 
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When  this  is  substituted  back  into  the  incremental  flow  rule  of  Eq, 


2.14,  the  expression  becomes  [16] 


^11*11  *11^22  ®11*33 


®22*11  *22*22  *22*33 


*33*11  *33*22  *33*33 
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where 


*11  “  ^11®11  *  '"l2®22 


*22  *  ’^12®11  °22 


*33  *  ^'’66°12 


In  symbolic  matrix  notation  this  is  written  as 


{deP}  «  [C]P  (do) 


(3.3) 


It  is  clear  that  the  incremental  plastic  flow  is  a  function  of 
the  current  stress  state,  effective  stress,  degree  of  work-hardening, 
and  the  current  increment  of  stresses  (do^j,  do22» 
interesting  to  note  that  the  incremental  plastic  strains  are  not  always 
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orthotropic  with  respect  to  the  incremental  stresses.  For  example, 
when  Oj2  non-zero,  it  is  possible  that  the  products  *11333 
®22®33  non-zero,  and  the  normal  incremental  strain  components  are 
then  coupled  to  the  incremental  shear  stress.  This  is  an  indication 
that  the  intuition  developed  from  elastic  analyses  may  not  always  aid 
in  understanding  orthotropic  plasticity. 

The  expression  for  the  total  incremental  stress-strain  relation¬ 
ship  is  obtained  by  summing  the  elastic  and  plastic  strain  increments 

{del  =  [c]  {do}  (3.4) 


where 

[c]  •  [c]'  ♦  [c]" 

The  [c]  matrix  relating  incremental  stresses  and  strains  will  be 
utilized  directly  when  establishing  the  element  stiffness  matrix. 

Numerical  Implementation  of  Loading  and  Unloading 

The  ANPLAST  program  is  designed  for  cyclic  loading  and  unloading. 
Therefore,  each  element  is  permitted  to  yield,  work-harden,  unload 
back  into  the  elastic  region,  and  then  re-load  back  to  the  new  yield 
surface.  This  transition  from  elastic  to  elastic-plastic  is  calculated 
on  the  element  level  and  is  not  constrained  by  the  direction  of 
loading  at  the  remote  location.  There  is  no  limit  to  the  number  of 
times  an  element  may  pass  through  the  elastic-plastic  transition.  This 
capability  will  make  possible  a  cycle  by  cycle  analysis  when 
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investigating  cyclic  fatigue.  At  the  present,  only  isotropic  hardening 
is  available,  but  other  hardening  theories  could  be  added  with  little 
difficulty. 

FORMULATION  OF  THE  CONSTANT  STRAIN  TRIANGLE 
The  linear  displacement  constant  strain  triangle  (CST)  [15,17]  was 
chosen  as  the  primary  finite  element  in  ANPLAST.  The  simplicity  of 
the  formulation  and  implementation  into  a  large  FORTRAN  program 
provides  considerable  justification  for  choosing  this  element. 

The  Strain-Nodal  Displacement  Relations 

A  typical  triangular  element  is  represented  in  Fig.  3.2,  along 
with  its  displaced  shape  [17].  Because  of  the  linear  displacement 
assumptions,  the  displacements  at  any  global  point  (x,y)  can  be 
expressed  as  the  following  linear  functions 

u(x,y)  =  u^  +  CjX  +  C2y 

(3.5) 

v(x,y)  =  +  C3X  +  c^y 

The  constants  c^,  Cg,  c^,  and  c^  can  be  obtained  in  terms  of  the  corner 
displacements  and  the  geometry  of  the  element. 

Cj,  Cg.  C3.  c^  -  F(u^.  V.,  Uj,  Vj,  u,^,  v^) 

The  three  global  strain  components  within  the  element  can  be  obtained 
from  the  assumed  displacement  field  as  follows 
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Figure  3.2  Dimensions  and  assumed  displacements  of  a  constant 
strain  triangle  finite  element. 
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^  =  c 
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M  =  C 
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(3.6) 


iH  + 

3y  3x 


C2  +  C3 


The  strains  may  be  represented  entirely  in  terms  of  nodal  displacements 
as  defined  in  Fig.  3.2. 
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where 


A  -  element  area 


or  in  symbolic  form 


{e}  »  [A]  {q} 


(3.7) 


where  {q}  are  the  nodal  displacements. 


Element  Stiffness  Matrix  Formulation 

The  incremental  element  stiffness  matrix  [K]  can  now  be  obtained 
by  use  of  Castigliano's  Theorem,  or  other  similar  energy  principles 
[18]  as 
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(3.8) 


[K]  =  II  [A]^[C][A]  dxdy 

Siilce  the  strains  and  stresses  are  constant  throughtout  the 
element,  all  three  of  the  matrices  may  be  brought  outside  the 
integrals,  and  the  area  integration  carried  out  explicitly.  Since  the 
element  is  triangular  and  has  straight  edges  the  area  integration  is 
direct  and  exact.  Note,  that  when  an  element  becomes  elastic-plastic, 
[C]  becomes  stress  dependent  and  consequently  so  does  [K].  This  is 
why  the  global  stiffness  matrix  needs  to  be  updated  at  the  beginning 
of  every  increment  when  plasticity  has  occurred. 

FORMULATION  OF  THE  ISOPARAMETRIC  ELEMENT 
It  is  generally  thought  that  the  isoparametric  element  offers 
some  advantage  in  efficiency  and  accuracy  compared  to  the  CST.  For 
this  reason,  a  quadratic  8-node  isoparametric  element  is  implemented 
into  ANPLAST.  The  descriptive  flow  diagram  and  subroutine 
descriptions  are  found  in  Figs.  3. 3-3. 4. 

The  Strain-Nodal  Displacement  Relations 

Isoparametric  elements  are  those  for  which  the  functional 
representation  of  deformational  behavior  is  employed  in  representation 
of  the  element  geometry  [18]. 


*  I  N^u^  y  »  I  N.v. 

u  *  I  N.x.  V  «  I  N.y. 


(3.9) 


where 
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Figure  3.3  Flow  Diagram  of  an  alternate  version  of  ANPLAST 
which  utilizes  an  isoparametric  quadralateral 
element. 
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INPUT: 

MESG: 

AMAT: 

JACOBS: 

LOAD: 

FORMK: 

ESMAT: 

TSTIFF: 

SSE: 

SSP: 

R0T4: 

ROT2 

INVERT: 

VCVTFG: 

LEQIPB: 

STRESS: 

OUTELM: 

OUTNOO: 

COEF: 


Figure 


Read  in  model  geometry,  boundary  conditions 
and  initialize  the  program  control  parameters. 

Provide  warning  and  error  messages. 

Form  the  element  strain-nodal  displacement  matrices. 

Form  Jacobian  matrix  and  determinant. 

Read  in  the  incremental  nodal  loads. 

Form  the  global  stiffness  matrix. 

Form  the  element  stiffness  matrix. 

Form  full  global  stiffness  matrix. 

Form  the  incremental  elastic  stress-strain  relations. 

Form  the  incremental  elastic-plastic  stress-strain 
relations. 

Angular  transformation  of  4th  order  tensor. 

Angular  transformation  of  2nd  order  tensor. 

Inversion  of  3x3  matrix. 

Transform  stiffness  matrix  to  band-symmetric  form. 
Solution  of  global  system  of  equations. 

Calculate  incremental  element  stresses  and  strains. 
Output  element  stresses  and  strains. 

Output  nodal  stresses  and  strains. 

Supply  coefficients  for  nodal  extrapolation. 


3.4  Subroutine  description  for  the  isoparametric 
element  version  of  ANPLAST. 
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x,y  =  geometry  of  element 
u,v  =  displacement  field 
=  shape  functions 

The  shape  functions  for  the  general  8-noded  isoparametric  element 
illustrated  in  Fig.  3.5  are 

-  J  (l-C)(l-n)(C+n+l)  N2  =  I  (l-c2)(l-n) 

N3  =  \  (l+e)(l-n)(C-n.l)  i  (lV)(l+C) 

(3.10) 

Ng  =  I  (l+C)(l+n)(C+n-l)  Ng  =  |  (l-c2)(l+n) 

N;  =  -  I  (l-C)(l+n)(C-n+l)  Ng  »  |  {i.n2)(i.c) 

The  strain-nodal  displacement  matrix  can  be  found  by  application 
of  the  basic  differential  definitions  of  strain  [18] 


or  in  symbolic  notation 

{e}  «  [A]  {q}  (3.11) 
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Figure  3.5  Local  coordinate  system  for  tfe  8-node  isoparametric 
quadralateral  element. 
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and  the  derivatives  are  calculated  by 
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where  [J]"^  is  the  inverse  Jacobian  transformation  matrix  which  can  be 
calculated  by 


[J]  = 


8  3N. 
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Element  Stiffness  Matrix  Formulation 

As  was  demonstrated  for  the  constant  strain  triangle  in  Eq.  3.8 
the  incremental  element  stiffness  matrix  [K]  can  be  obtained  by  use  of 
Castigliano's  Theorem 


tK]  .  fj 


[A]  [C][A]  dxdy 


(3.12) 


Inspection  of  Eq.  3.11  reveals  the  strains  and  stresses  are  not 
constant  throughout  the  element.  Therefore,  the  matrix  product 
[A]^[C][A]  cannot  be  taken  outside  the  integral.  Consequently,  a 
numerical  integ'*ation  scheme  must  be  Implemented  in  the  program. 
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Numerical  Integration  of  Element  Matrices 
Using  the  following  transformation. 


dxdy  *  det  [J]  dcdn 


the  element  stiffness  matrix  now  becomes 

1  1 


[K]  -  I  I  [A]^CC][A]  |J1  dcdn 
-1  -1 


(3.13) 


Now,  by  using  the  Gaussian  quadrature  technique,  the  element 
stiffness  matrix  can  be  calculated  by  summing  the  matrix  quantity 
CA]^[C][A]  at  the  four  selected  Gauss  points  (c  ■  ±.577350, 
n  *  t. 577350)  and  multiplying  each  by  the  appropriate  weighting  factor 
(which  in  this  case  is  unity).  In  suwnary,  the  stiffness  matrix  can 
be  represented  in  symbolic  form  as 

[K]  ■  I  [A]^[C][A]  iJ|  (3.14) 

i-1  ’ 

w.j  =  1.0  (weighting  function) 

It  should  be  noted  that  for  non>linear  plasticity  analyses, 
partial  yielding  may  occur  in  a  given  element.  One  particular  Gauss 
point  may  be  yielded  while  the  others  are  still  elastic.  This  will 
not  present  a  problem  when  formulating  the  element  stiffness  matrix 
since  the  elastic-plastic  constitutive  relations  [C]^  can  be  utilized 
when  yielding  occurs  at  a  particular  Gauss  point,  and  then  the  element 
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stiffness  matrix  summed  in  the  usual  way. 


PRE-PROCESSING  AND  POST-PROCESSING  FEATURES 
To  make  the  utilization  of  ANPLAST  more  efficient  and  convenient 
a  number  of  pre-  and  post-processing  software  have  been  developed. 

Mesh  Generation  and  Plotting  Features 

The  creation  of  a  geometrical  mesh  is  the  first  step  in  a  finite 
element  analysis  and  can  often  be  the  most  time  consuming.  Therefore, 
four  generic  mesh  generators  have  been  developed  for  rectangular, 
tapered  quadrilateral,  radial  fan,  and  circular  hole  geometries.  An 
example  of  each  type  of  mesh  is  illustrated  in  Fig.  3.6. 

All  of  the  programs  operate  on  an  interactive  basis,  prompting  the 
user  for  appropriate  input  parameters  such  as  maximum  and  minimum 
dimensions  as  well  as  mesh  density.  The  mesh  generators  create  two 
files,  one  containing  the  nodal  coordinates  and  the  other  the  element 
connectivity.  The  starting  node  and  element  numbers  are  required  so 
that  successive  files  can  be  combined  for  constructing  a  more  complex 
geometry . 

Bandwidth  Optimization  by  Nodal  Sweepinq 

A  nodal  sweeping  technique  was  developed  to  join  together  the 
various  composite  meshes  into  one  integral  mesh.  The  interactive  nodal 
sweeping  program  works  by  passing  a  plane  or  sphere  through  the 
composite  mesh  and  fusing  it  into  one  by  removing  any  duplicate  modes. 
This  sweeping  process  also  renumbers  the  nodal  connectivity  list  and 
alters  the  stiffness  matrix  bandwidth.  A  judicious  chose  of  sweeping 
orientation  can  lead  to  a  significant  reduction  in  bandwidth. 


BIASED  QUADRILATERAL 


TAPERED  QUADRILATERAL 
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WS¥ji¥M 

WSSMWM 


«« 


RADIAL  FAN  CIRCULAR  HOLE 

(CRACK  TIP) 


Figure  3.6  A  collage  of  finite  element  meshes  created 
by  the  ANPLAST  pre- processors.  Each  mesh 
design  demonstrates  biasing  of  element 
density. 
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Considerable  gains  in  analysis  efficiency  have  been  evident  since 
developing  the  nodal  sweeping  program. 

Nodal  Averaging  of  Stresses  and  Strains 

Solution  of  the  finite  element  equations  provides  element  stresses 
and  strains.  For  the  case  of  the  isoparametric  element,  the  stresses 
are  provided  at  the  Gaussian  integration  points.  The  constant  strain 
triangle  assumes  the  stress  to  be  constant  across  the  area  of  the 
element.  Unfortunately,  the  element  or  Gauss  point  locations  usually 
have  no  special  geometrical  significance,  while  nodes  often  fall  on 
convenient  lines  through  the  mesh.  Therefore,  nodal  averaging 
techniques  were  developed. 

For  both  the  triangular  and  quadrilateral  elements,  the  nodal 
values  of  stresses  are  found  by  averaging  the  contribution  of  stress 
from  each  connecting  element  [18].  For  the  case  of  the  isoparametric 
element,  individual  element  nodal  stresses  are  calculated  by  extra¬ 
polation  of  Gauss  point  stresses  using  a  modified  form  of  the  element 
shape  functions  [19].  The  nodal  stresses  can  be  provided  directly  by 
ANPLAST  in  either  Cartesian  or  polar  coordinate  systems. 

Presentation  of  Results  by  Contour  Plotting 

The  vast  magnitude  of  results  generated  by  a  finite  element 
analysis  makes  necessary  innovati-ve  techniques  for  post  processing. 

One  such  technique  is  the  contouring  of  stresses,  strains,  displace¬ 
ments,  or  any  combination  of  these  quantities.  In  this  particular 
case,  the  contouring  is  based  on  interpolation  of  nodal  quantities, 
and  plotted  using  the  NCAR  package  available  on  the  Purdue  University 
Computing  Center  system. 
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Sections  4,  5,  and  6  will  provide  examples  of  the  use  of  the  post¬ 
processing  programs  developed  for  this  research. 

NUMERICAL  EXTENSION  OF  LAMINATION  THEORY 
LAMINATE  Computer  Program 

The  macroscopic  laminate  behavior  described  in  Section  2  is 
adapted  for  use  in  an  interactive  FORTRAN  77  computer  program  called 
LAMINATE.  The  purpose  of  the  LAMINATE  program  is  to  study  the  multi - 
axial  response  of  elastic-plastic  laminates.  The  incremental  laminate 
load  deformation  equations  can  be  written  as 


where 


f  0^ 

dN 

de° 

X 

X 

>  •  [Nj]  < 

"'J 

dN 

dv^,. 

L  "V 

(3.15) 


LAMINATE  logic  has  been  designed  to  parallel  that  of  ANPLAST  so 

that  loading  increments  (dN  ,  dN  ,  dN  )  are  applied  and  the  resulting 

X  y  xy 

mid-plane  strains  calculated.  This  is  accomplished  by  inverting  [A^j] 
after  it  is  numerically  evaluated. 

As  the  incremental  analysis  progresses,  the  individual  layers  can 
become  elastic-plastic,  and  their  contribution  ente-s  through  the 
formulation  of 

I J  ^ 
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The  LAMINATE  program  can  analyze  a  laminate  consisting  up  to  10 
orthotropic  elastic-plastic  layers.  Each  layer  is  characterized 
individually,  and  obeys  its  own  yield  and  work-hardening  law.  The 
program  output  consists  of  individual  layer  stresses  and  overall 
strains. 

Modifications  to  ANPLAST 

The  ANPLAST  finite  element  program  was  modified  to  allow  for  ortho¬ 
tropic  elastic-plastic  laminated  material  response.  Basically,  the 
lamination  theory  developed  for  the  LAMINATE  program  is  applied  to 
ANPLAST.  The  primary  modifications  occur  in  the  MAIN  program  and  in 
subroutine  ESMAT  where  it  is  necessary  to  store  each  of  the  individual 
layer  stresses  and  check  for  yielding  in  each  layer  of  each  element. 

When  formulating  the  element  stiffness  matrix  it  is  necessary  to  sum 
over  each  layer  through  the  thickness  of  the  laminate.  Only  minor 
modifications  are  made  to  the  standard  ANPLAST  input  structure.  It  is 
interesting  to  note  that  analyzing  a  single  layer  is  a  special  sub-set 
of  lamination  theory,  and  therefore  this  version  of  ANPLAST  can  be  used 
for  all  analyses. 


COMPUTER  PROGRAM  PORTABILITY 

All  of  the  computational  programs  developed  in  this  section  have 
been  designed  to  be  portable.  None  of  the  programs  utilize  specialized 
utility  routines,  such  as  those  from  IMSL.  ANPLAST  and  LAMINATE  can 
be  easily  executed  on  interactive  as  well  as  batch  processors.  At 
Purdue,  both  of  these  programs  have  been  successfully  run  on  the  VAX 
11/780,  CDC6600,  and  CYBER  205.  The  majority  of  pre-processing  and 


post-processing  programs  are  best  run  on  an  interactive  system  as  this 
provides  the  greatest  convenience  for  the  user. 

Since  all  of  the  programs  utilize  only  standard  FORTRAN  77  library 
routines,  they  are  completely  portable  to  other  systems.  The  only 
limitation  will  be  the  amount  of  memory  made  available  to  the  programs. 
As  the  problems  become  larger,  computers  with  a  virtual  memory  system 
appear  to  offer  the  only  alternative  unless  out  of  core  solution 
algorithms  are  developed. 


f 
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SECTION  4  ELASTOPLASTIC  STUDY  OF  A  SHEET  WITH  A  CIRCULAR  HOLE 


INTRODUCTION 

A  thin  sheet  with  a  circular  hole  is  analyzed,  because  it  is  one 
of  the  standard  benchmark  problems.  It  is  expected  that  many  of  the 
trends  discovered  for  the  hole  will  £fpply  to  more  severe  discontin¬ 
uities,  such  as  cracked  panels.  Additionally,  there  is  a  body  of 
literature  available  containing  some  theoretical  and  experimental 
elastoplastic  solutions. 

The  ANPLAST  program  will  be  used  to  investigate  the  isotropic 
elastoplastic  response  of  a  sheet  with  a  circular  hole  and  comparisons 
made  with  some  of  the  solutions  in  the  literature.  This  study  will 
then  be  extended  to  orthotropic  plastic  flow,  by  conducting  a  parametric 
study  of  the  plasticity  coefficients.  The  numerical  results  are 
compared  to  an  elastoplastic  experimental  study  of  a  metal  matrix  fiber 
reinforced  sheet.  A  problem  with  a  remote  complex  loading  path  is  also 
considered. 


SHEET  WITH  ISOTROPIC  MATERIAL  PROPERTIES 
Elastic  Loading:  Comparison  with  Theoretical  Solution 

A  thin  sheet  with  isotropic  elastic  material  properties  was  modeled 
by  ANPLAST  and  analyzed  for  the  case  of  a  uniaxial ly  applied  remote 
load.  The  finite  element  mesh  used  is  shown  in  Fig.  4.1.  The  fixed 
grip  end  condition  simulated  by  inputting  very  stiff  material 
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Figure  4.1  Fir ite  element  mesh  of  a  sheet  with  a  circular 
hole  subjected  to  a  uniaxial  remote  load  with 
fixed  grip  end  conditions. 
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properties  to  the  top  row  of  elements,  and  restraining  rotation  by 
applying  two  additional  boundary  restrains. 

In  the  following  discussions  notational  definitions  are  as  shown 
in  Fig.  ^,2.  A  comparison  of  the  nodal  stresses  and  theoretical  solution 
along  a  horizontal  radial  line  (e  =  0°)  is  presented  in  Fig.  4.3.  The 
theoretical  solution  obtained  from  the  computer  program  ELLPLAS 
[12]  was  also  used  for  the  orthotropic  analysis.  It  is  clear 
that  the  numerical  model  agrees  quite  well  with  theory.  The  small 
discrepancies  are  attributed  to  the  theoretical  solution  failing  to 
account  for  the  finite  dimensions  of  the  sheet. 

Figure  4.4  illustrates  the  elastic  effective  stress  (a)  distribu¬ 
tion  along  three  radial  lines  (e  =  O'^,  45°,  90°),  and  again  is  in  close 
agreement  with  theory.  The  effective  stress  is  of  special  interest 
since  it  indicates  the  onset  of  plasticity.  Figure  4.4  clearly 
suggests  a  strong  tendency  of  yielding  to  begin  along  the  horizontal 
symmetry  line  (9  =  0°)  at  the  edge  of  the  hole.  However,  at  about  1.5 
radii  from  the  hole,  o  along  the  45°  orientation  begins  to  dominate. 
This  suggests  that  growth  of  any  potential  plasticity  may  begin 
horizontally,  but  will  eventually  produce  a  lobe  in  the  45°  direction. 
Of  course,  the  results  presented  are  for  a  purely  elastic  analysis  and 
do  not  account  for  redistribution  of  stresses  during  plastic  flow. 

An  iso-stress  contour  of  the  effective  stress  is  presented  in 
Fig.  4.5.  This  unique  contour  pictorially  represents  the  elastic 
prediction  of  the  plastic  zone  shape  at  the  onset  of  yielding.  Since 
the  contours  are  labeled  with  the  appropriate  values  of  o,  this  figure 
also  provides  the  location  of  initial  yield,  and  the  direction  of 


0^  =  11  Material  Direction 
0  =  Angular  Orientation  of  Element 
r  =  Radial  Distance 


Figure  4.2  Local  and  global  coordinate  systems  for  a 
sheet  with  a  circular  hole. 
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3  Comparison  of  ANPLAST  nodal  stresses  with  the 
ELLPLAS  theoretical  solution  for  a  sheet  with 
a  circular  hole  and  isotropic  material 
properties. 
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Figure  4.4  C<»par1son  of  ANPLAST  effective  stress  distribution 
with  the  ELLPLAS  theoretical  solution  for  an 
isotropic  sheet  with  a  circular  hole. 
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5  Effective  stress  contour  for  the  sheet  with  isotropic 
material  properties  subjected  to  elastic  loading. 
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subsequent  plastic  flow.  The  contour  plot  is  consistent  with  the 
previous  discussion  involving  the  effective  stress  distribution. 

ANPLAST  Replication  of  Ratnberq-Osqood  Stress-Strain  Behavior 

The  isotropic  elastic-plastic  constitutive  model  chosen  was  that 
of  a  Ramberg-Osgood  work-hardening  type  material,  whose  uniaxial 
behavior  is  represented  as  [20] 

e  •  f  [1  *1  (f-)"  ]  (4.1) 

where 

E  -  elastic  modulus 

Oy  »  uniaxial  yield  stress 

n  »  work-hardening  power 

0  *  uniaxial  stress 

The  elastic  and  plastic  strain  components  can  be  separated,  and 
can  be  written. 


I 


(^.2) 

V 

i 


For  a  multiaxial  stress  state,  the  uniaxial  stress  quantities  are 
replaced  by  the  effective  values,  as  defined  in  Section  2.  The 
effective  plastic  strain  can  now  be  written  as  a  function  of  the 
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effective  stress 


eP  *  I  f  (I-)"  ^  (4.3) 

The  ANPLAST  program  is  designed  so  that  any  general  work-hardening 
law  can  be  utilized.  Therefore,  according  to  Eq.  2.14,  the 
slope  of  the  work-hardening  relationship  must  be  explicity  defined. 

For  the  Ramberg-Osgood  work-hardening  formulation,  this  is  easily 
accomplished  by  differentiating  Eq.  4.3  and  inverting,  to  give 


deP 


—  (S-) 
en 


-  (l-n) 


(4.4) 


where 

B  =  curve  fitting  parameter 
o  =  effective  stress 
Oy  *  effective  yield  stress 

In  addition,  the  plasticity  coefficients  are  assumed  to  have  the 
following  values  for  an  isotropic  material 


r^2  ‘  -0.50 
rgg  *  1.50 
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ANPLAST  results  for  a  uniaxial 1y  loaded  isotropic  elastic-plastic 
specimen  are  plotted  in  Fig.  4.6  and  compared  to  the  ideal  Ramberg- 
Osgood  material  response  represented  by  Eq.  4.1.  The  small  deviations 
in  the  comparison  are  due  to  the  nature  of  the  elastic-plastic  finite 
element  analysis.  While  Eq.  4.1  represents  a  smooth  continuous 
behavior,  ANPLAST  assumes  elastic,  then  elastic-plastic  response 
characterized  by  a  distinct  yield  point.  These  results  show  that  the 
Ramberg-Osgood  form  of  the  stress-strain  behavior  is  adequately 
replicated.  This  is  important  for  the  later  comparisons. 

Elastoplastic  Response 

Isotropic  elastic-plastic  response  of  the  circular  hole  problem 
is  of  interest  for  two  important  reasons.  First,  it  is  one  of  the  few 
non-trivial  problems  where  adequate  theoretical  and  experimental 
solutions  exist.  Second,  it  provides  a  benchmark  from  which  ortho¬ 
tropic  eKstoplastic  behavior  can  be  gauged.  The  isotropic  elastic- 
plastic  constitutive  model  chosen  is  that  of  a  Ramberg-Osgood  work¬ 
hardening  type  material,  as  discussed  in  the  previous  section. 

An  incremental  elastic-plastic  finite  element  analysis  of  the 
circular  hole  was  conducted.  The  analysis  was  conducted  by  first 
loading  the  model  elastically  until  the  first  element  yielded,  then 
taking  10  equal  remote  load  steps  of  3000  ps i .  As  predicted  by  the 
elastic  analysis,  yielding  occurred  at  the  hole  along  the  horizontal 
symmetry  line.  Yield  of  the  first  few  elements  took  place  at  a 
remote  load  of  17000  psi.  This  is  reasonable  since  the  elastic  stress 
concentration  factor  is  approximately  3  and  the  uniaxial  yield  stress 
is  50000  psi . 
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A  visualization  of  plastic  zone  growth  is  illustrated  in  Fig.  4.7. 
This  depiction  of  plastic  zone  growth  is  accomplished  by  placing  an 
asterisk  (*)  at  the  centroid  of  each  element  which  is  currently 
yielded.  The  plastic  zone  initially  progresses  horizontally  then  forms 
a  lobe  at  around  e  =  45°  at  the  higher  loads.  It  is  interesting  to 
note  that  net  section  yielding  does  not  occur  at  the  horizontal 
symmetry  line,  but  rather  the  lobe  at  e  *  45°  extends  across  the 
specimen. 

An  experimental  analysis  of  the  plastic  zone  growth  in  an 
isotropic  elastic-plastic  material  was  conducted  by  P.S.  Theocaris  and 
E.  Marietos  in  1964  [21].  Utilizing  the  methods  of  birefringent 
coatings  and  electrical  analogy,  they  determined  the  elastoplastic 
response  of  an  aluminum  alloy  (57S)  under  remote  uniaxial  loading.  The 
plastic  zone  boundaries  which  they  found  during  subsequent  load 
increments  (I-VIII),  are  shown  in  Fig.  4.8.  The  growth  of  these 
plastic  zones  progress  in  the  same  fashion  as  that  calculated  by 
ANPLAST  in  the  preceding  analysis.  Though  a  direct  quantitative 
comparison  cannot  be  offered  due  to  the  different  material  properties, 
geometry,  and  loading,  the  character  of  the  results  are  the  same. 

Budiansky  and  Vidensek  [20]  arrived  at  an  approximate  theoretical 
solution  for  the  problem  of  a  circular  hole  in  an  infinite  plate  for  an 
isotropic  elastoplastic  material.  They  also  assumed  the  elastic- 
plastic  work-hardening  was  of  the  Ramberg-Osgood  type.  Using  this 
solution,  the  stresses  and  strains  can  be  calculated  at  any  point  in 
the  body. 
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Figure  4.7 


Plastic  zone  growth  for  an  isotropic  sheet 
subjected  to  elastic-plastic  loading. 
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Figure  4.8  Plastic  zone  growth  experimentally  detained  by  Theocaris  and  Marietos  [21] 
using  a  birefringent  coating  on  an  isotropic  sheet,  then  subjected  to 
elastic-plastic  loading. 


The  ANPLAST  nodal  stresses  were  plotted  along  three  radial  lines 
(9  =  0°,  45°,  90°)  and  compared  to  the  Budiansky-Vidensek  solution  in 
Fig.  4.9  for  two  different  remote  loads.  The  effective  stress  distri¬ 
bution  for  the  same  two  load  levels  are  shown  in  Fig.  4. 10.  The 
finite  element  results  are  in  good  agreement  with  the  approximate 
theoretical  solution,  with  the  greatest  deviation  occurring  at  the 
higher  remote  load.  Part  of  the  discrepancy  can  be  related  to  the 
theoretical  solution  failing  to  account  for  the  finiteness  of  the 
specimen  modeled  by  ANPLAST. 

In  1958,  Frocht  and  Thomson  [22]  studied  the  isotropic  elasto- 
plastic  response  of  a  sheet  with  a  circular  hole  using  the  method  of 
photoplasticity.  The  results  for  a  Ramberg-Osgood  type  work-hardening 
material  are  presented  in  Fig.  4.\1  alongside  the  ANPLAST  and 
Budiansky-Vidensek  solutions.  In  general,  the  photoplasticity  results 
fall  onto  the  same  distribution  as  the  other  two  solutions.  They  also 
fall  above  the  Budiansky-Vidensek  solution  at  a  radius  of  about  1.5  - 
2.0  inches,  as  did  the  finite  element  results.  It  should  be  noted  that 
the  actual  dimensions  of  Thompson's  test  specimen  were  slightly 
different  than  that  modeled  by  ANPLAST.  The  results  consequently  were 
normalized  with  respect  to  the  hole  radius  and  remote  load. 

The  isotropic  elastoplastic  analysis  capability  of  the  ANPLAST 
program  has  been  verified. 

Comparison  of  the  Isoparametric  Quadralateral  and  Constant  Strain 
Triangle  Elements 

It  has  been  suggested  by  a  number  of  authors  [19,23  that  the 
higher  order  isoparametric  elements  may  provide  superior  performance 
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Figure  4.9  Comparison  of  ANPLAST  nodal  stresses  with  the 
Budiansky-Vidensek  [20]  theoretical  solution 
for  an  isotropic  sheet  with  a  circular  hole 
subjected  to  elastic-plastic  loading. 
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Figure  4.10  Comparison  of  ANPLAST  effective  stresses  with  the 
Budiansky-Vidensek  [20]  solution  for  an  isotropic 
sheet  with  a  circular  hole  subjected  to  elastic- 
plastic  loading. 
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Figure  4.11  Comparison  of  ANPLAST  stresses  with  an  experimental 
photoplasticity  solution  by  Frocht  and  Thomson  [22] 
for  the  isotropic  sheet  with  a  circular  hole 
subjected  to  elastic-plastic  loading. 
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for  certain  plasticity  analyses.  The  8-noded  quadralateral  with  a 
4-point  Gaussian  integration  scheme  developed  in  Section  3  was  used  to 
analyze  the  elastoplastic  response  of  the  isotropic  sheet  with  a 
circular  hole.  The  plastic  zone  visualization  for  this  particular 
element  is  handled  by  placing  an  asterisk  (*)  at  any  Gauss  point  which 
has  exceeded  the  yield  stress. 

The  mesh  used  in  this  analysis*  shown  in  Fig.  .12,  is  of  the 

same  geometry  as  in  the  previous  analysis  conducted  with  the  CST.  In 
addition,  the  material  properties  are  also  identical.  The  elastic 
stresses  are  plotted  in  Fig.  4.13  along  the  horizontal  synmetry  line 
(e  =  0°).  Though  there  are  more  degrees  of  freedom  for  this  mesh  than 
the  corresponding  CST  mesh,  the  results  are  no  better.  In  fact,  the 
nodal  stress  (o^y)  on  the  hole  boundary  displays  considerably  more 
deviation.  This  is  probably  due  to  the  particular  extrapolation 
technique  [19]  used  to  derive  nodal  stresses  from  the  Gauss  points. 

The  plastic  zone  progression  is  illustrated  in  Fig.  4.14  for  the 
same  remote  stress  levels  as  was  presented  for  the  CST.  The  results 
for  the  quadrilateral  element  demonstrate  excellent  correlation  with 
the  CST  to  the  point  where  any  discrepancies  are  almost  indistinguish¬ 
able.  These  findings  suggest  that  there  is  no  practical  advantage  for 
selecting  the  isoparametric  element  over  the  CST,  at  least  as  far  as 
isotropic  elastoplastic  analyses  are  concerned.  Therefore,  due  to  the 
simplicity  and  efficiency  of  the  CST,  this  element  will  be  the  primary 
analysis  tool  for  this  research. 
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Figure  4.13  Isoparametric  quadralateral  stress  distribution 
for  an  isotropic  sheet  with  a  circular  hole 
subjected  to  elastic  loading. 
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Figure  4.14  Comparison  of  plastic  zone  growth  for  the  isoparametric 
quadra! ateral  and  CST  elements  for  an  isotropic  sheet 
subjected  to  elastic-plastic  loading. 


86 


SHEET  WITH  ORTHOTROPIC  MATERIAL  PROPERTIES 
Elastic  Loading:  Comparison  with  Theoretical  Solution 

An  elastic  analysis  of  the  sheet  with  a  circular  hole  was  conducted 
for  two  orthotropic  materials.  The  first,  Boron-Aluminum,  demonstrates 
mild  orthotropy,  while  the  latter,  Graphite-Epoxy,  exhibits  rather 
strong  anisotropy.  Two  analyses  were  conducted  for  each  material,  one 
with  the  stiff  material  direction  being  oriented  horizontally 
(e^  =  0°),  the  other  with  it  rotated  90°  (e^  *  90°). 

The  nodal  stresses  along  the  horizontal  symmetry  line  are  plotted 
against  the  infinite  sheet  theoretical  solution  generated  by  the 
program  ELLPLAS  in  Figs.  4.15-4.16.  The  B-Al  results  are  similar  to 
the  Isotropic  material  which  is  what  one  would  expect  since  this 
material  is  only  weakly  orthotropic.  The  6-Ep  findings  highlight  the 
effect  of  strong  anisotropy  on  the  stress  distributions.  When  the 
anisotropy  is  strong  and  the  stiff  material  direction  is  vertical,  the 
gradient  of  in  the  vicinity  of  the  hole  becomes  quite  steep.  Also, 
this  orientation  appears  to  significantly  increase  the  elastic  stress 
concentration  factor,  while  the  opposite  trend  is  true  for  »  0°. 

The  finite  element  results  are  in  good  agreement  with  the 
theoretical  solution  in  most  locations.  For  the  G-Ep(e^  *  0°) 
specimen,  a  considerable  deviation  is  apparent  for  away  from  the 
vicinity  of  the  hole.  This  can  be  accounted  for  by  recalling  that  the 
theoretical  solution  is  for  the  infinite  sheet.  While  the  finite  size 
effect  did  not  manifest  itself  noticeably  for  the  isotropic  material, 
it  appears  that  strong  anisotropy  can  increase  the  finitehess  of  a 
particular  specimen.  It  is  clear  from  these  results  that  the 
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rigure  4.15  Comparison  of  the  ANPLAST  stress  distribution  with 
the  ELLPLAS  theoretical  solution  for  a  mildly 
orthotropic  (B-Al)  sheet  with  a  circular  hole 
subjected  to  elastic  loading. 
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gure  ‘».16  Comparison  of  the  ANPLAST  stress  distribution  with  the  ELLPLAS  theoretical  sol 
for  a  strongly  orthotropic  (G-Ep)  sheet  with  a  circular  hole  subjected  to  elas 
loading. 


orthotropic  elastic  capabilities  of  ANPLAST  have  been  verified. 

A  Study  of  the  Plasticity  Coefficients  Effect  on  Plastic  Zone  Growth 
It  is  the  intention  of  this  study  to  investigate  the  effect  of 
introducing  anisotropy  into  the  plastic  flow  relations  for  the  sheet 
with  a  circular  hole.  Anisotropy  is  brought  into  the  problem  through 
the  plasticity  coefficients,  r^j,  r^g.  The  previous  elastoplastic 
analysis  assumed  isotropic  values  of  the  plasticity  coefficients.  In 
this  investigation,  the  plasticity  coefficients  r^^^,  and  rgg  are  varied 
to  determine  how  they  affect  yielding  and  subsequent  plastic  zone 
progression.  Since  the  anisotropy  of  the  plasticity  is  of  primary 
interest,  the  examples  in  this  section  assume  elastic  isotropy. 

The  first  analysis  involves  the  variation  of  the  plasticity 
coefficient  r^  from  its  isotropic  value  of  unity.  The  assumed  value 
of  r^  is  chosen  to  be  0.20.  The  uniaxial  yield  stress  in  the  strong 
direction  for  this  material  is  2.23  times  greater  than  in  the  weak 
direction.  At  the  same  time,  it  is  necessary  to  reduce  the  magnitude 
or  rj2  -0.50  (isotropic)  to  -0.20,  to  prevent  the  anisotropic 
effective  stress  from  becoming  negative.  This  finding  suggests  that 
there  are  additional  restrictions  on  the  plasticity  coefficients  than 
were  mentioned  in  Reference  [12],  since  a  negative  a  would  indicate 
negative  incremental  work,  which  is  physically  unreasonable. 

When  rjj  is  not  unity,  planar  anisotropy  is  introduced  in  the 
sheet,  resulting  in  directional  sensitivity  of  the  plastic  flow.  Two 
analyses  were  conducted,  one  with  the  strong  direction  oriented 
horizontally  (e^  =  0°)  and  the  other  with  the  weak  direction  horizon¬ 
tally  (e^  *  90°).  THe  Ramberg-Osgood  work-hardening  law  and  all  the 
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associated  coefficients  are  as  for  the  isotropic  elastoplastic 
analysis. 

The  plastic  zone  visualization  for  both  analyses  is  pictured  in 
Figs.  4.17-4.18.  In  comparison  to  the  isotropic  case,  it  is  clear  that 
the  coefficient  r^^  plays  a  dominant  role  in  plastic  zone  progression. 
While  the  plots  for  =  0°  are  somewhat  similar  to  those  for  the  isotropi 
material,  this  is  clearly  not  the  case  for  *  90°.  When  the  strong 
direction  is  oriented  parallel  to  the  load,  it  is  evident  that  the 
plastic  zone  progresses  in  the  vertical  direction,  and  does  not 
demonstrate  the  tendency  to  produce  net  section  yielding.  This 
behavior  is  understandable,  since  the  major  stress  component  (Oyy)  is 
parallel  to  the  strong  principal  material  direction,  and  is  not  as 
likely  to  cause  yielding. 

Yielding  at  0  *  90°  is  apparent  for  *  90°,  yet  is  absent  for 
9^*  0°.  Even  the  isotropic  example  eventually  exhibited  yielding  at 
this  location.  The  explanation  for  this  behavior  is  directly  related 
to  the  plastic  anisotropy. 

An  iso-stress  contour  of  effective  stress  at  the  onset  of  yield 
is  presented  for  both  analyses  in  Fig.  4.19.  These  plots  show  that 
the  initial  progression  of  the  plastic  zone  in  the  vicinity  of  the 
hole  is  strongly  correlated  with  the  stress  distribution  for  the 
elastic  solution.  These  contours  also  highlight  the  strong  gradient 
of  o  at  the  horizontal  hole  perimeter  (e  «  0°). 

Effective  stress  contours  for  =  0°,  and  ®  90°  at  remote 
loads  of  34000  psi  and  76000  psi,  respectively  are  shown  in  Fig.  4.20. 

Both  of  these  load  levels  represent  a  state  where  considerable 
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Figure  4.17  Plastic  zone  growth  for  an  orthotropic  sheet  subjected 
to  elastic-plastic  loading.  The  strong  material 
direction  is  oriented  perpendicular  to  the  remote 
load. 
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Figure  4.18  Plastic  zone  growth  for  an  orthotropic  sheet  subjected 
to  elastic-plastic  loading.  The  strong  material 
direction  Is  oriented  parallel  to  the  remote  load. 
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Figure  4.19  Effective  stress  contours  for  three  orthotropic  sheets 
at  the  onset  of  yielding. 


20  Effective  stress  contours  for  three  orthotropic  sheets  after  consi 


plasticity  has  occurred.  These  contours  again  accurately  portray  the 
outline  of  the  elastic  elastic-plastic  boundary,  and  are  strikingly 
similar  to  the  elastic  contours.  One  could  hypothesize  that  the  shape 
of  these  contours  will  not  deviate  to  any  great  extent,  until  the 
plasticity  has  grown  to  such  a  scale  that  it  interacts  with  the  finite 
boundaries. 

That  plastic  work-hardening  has  occurred  behind  the  initial  yield 
boundary  can  be  detected  by  noticing  an  increase  of  some  contours  above 
the  initial  yield  (a^).  Although  this  study  has  not  approached  the 
subject  of  work-hardening  and  its  effect  on  plastic  zone  growth,  it 
can  be  assumed  that  this  parameter  will  control  the  rate  at  which  the 
growth  occurs.  It  is  likely  that  the  degree  of  work-hardening  will  not 
affect  the  shape  of  the  plastic  zone,  rather  that  role  is  reserved  for 
the  plasticity  coefficients. 

A  similar  investigation  of  the  plasticity  shear  coefficient  rgg 
wasconducted  by  setting  r^^  back  to  its  isotropic  value  (1.0)  and 
letting  stay  at  -.20.  The  coefficient  rgg  was  doubled  from  its 
previous  value  and  assumed  to  be  3.0.  Both  *  0°  and  9^  ■  90°  were 
studied,  but  the  analyses  indicated  that  these  two  cases  produce  nearly 
identical  results.  This  finding  makes  sense  if  one  considers  that 
rgg  multiplies  only  the  shear  term,  and  the  transformation  of  the 
stresses  to  either  material  coordinate  system  orientation  will  produce 
a  shear  of  the  same  absolute  magnitude.  Therefore,  only  one  analysis 
will  be  discussed. 

The  plastic  zone  visualization  for  this  analysis  is  presented  in 
Fig.  4.21  and  the  results  are  almost  identical  to  those  of  Fig.  4.17. 
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Figure  4.21  Plastic  zone  growth  for  an  orthotropic  sheet 
(rgg  ■  3.0)  under  elastic-plastic  loading. 

The  orientation  of  the  principal  material 
directions  does  not  affect  this  illustration. 
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One  of  the  subtle  differences  is  a  slight  tendency  for  the  plastic 
zone  to  drift  a  bit  more  in  the  +6  direction  at  the  lower  remote  loads. 
This  can  be  associated  with  the  accentuation  of  the  shear  stress  which 
is  a  maximum  at  approximately  e  =  45°. 

The  effective  stress  contour  at  the  inception  of  yield  is 
presented  in  Fig.  4.19.  This  contour  is  nearly  identical  to  the 
previous  study  of  r^^  for  =  0°.  In  summary,  it  appears  that  the 
effect  of  rgg  on  plastic  zone  growth  for  this  particular  geometry, 
loading,  and  principal  material  orientation  is  negligible. 

A  Study  of  Unloading  and  the  Resulting  Residual  Stresses 

A  most  interesting  aspect  of  studying  plasticity  problems  is  the 
effect  of  load  reversal  or  unloading  and  the  resulting  re"’'dual 
stresses  which  are  created.  When  a  sheet  with  a  circular  hole  is 
loaded  beyond  its  elastic  limit  a  plastic  zone  results.  The  material 
inside  this  zone  is  permanently  deformed,  and  upon  load  reversal  the 
elastic  material  surrounding  it  tends  to  clamp  down  and  produce 
significant  compressive  stresses.  It  is  the  purpose  of  this  section 
to  determine  if  the  presence  of  plastic  anisotropy  will  have  a 
prominent  effect  on  the  residual  stress  distributions. 

The  first  unloading  case  Involves  the  Isotropic  elastoplastic 
material  analyzed  earlier  in  this  section,  and  will  serve  as  the  basis 
for  comparison.  In  all  cases,  the  maximum  remote  load  was  chosen  to 
be  two  times  the  load  required  to  initiate  yielding.  The  result  of 
this  was  a  maximum  remote  load  of  34000  psi  for  all  examples  except 
for  the  anisotropic  material  where  r^j  «0.20  and  ■  90°. 
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The  unloading  was  carried  out  by  first  taking  a  small  unloading  step 
which  resulted  in  all  elements  becoming  elastic  (falling  inside  their 
yield  surfaces).  With  this  accomplished,  and  the  global  stiffness 
matrix  reassembled,  one  large  unloading  step  to  zero  was  taken.  In  no 
case  did  any  elements  yield  again.  If  yielding  had  occurred, 
incremental  unloading  would  have  been  required. 

The  residual  stress  distribution  for  each  analysis  is  presented 
in  Figs.  4.22.  Since  the  maximum  load  and  plastic  zone  size  are  not 
the  same  for  each  case,  a  direct  comparison  of  the  stress  levels  is 
not  appropriate.  But  it  is  clear  that  the  character  of  the  results 
are  the  same.  The  vertical  stress  component,  demonstrates  a 
significant  negative  residual  stress  component  in  the  vicinity  of  the 
hole  perimeter  while  the  other  stresses  are  negligible.  In  each  case 
a  compressive  residual  stress  (o  )  is  present  out  to  approximately 
i"  to  1"  away  from  the  hole.  In  general,  this  brief  study  suggests 
that  planar  anisotropy  does  not  have  a  significant  effect  on  the 
residual  stresses  due  to  unloading. 

A  COMPLEX  LOADING  PATH 

All  of  the  previous  analyses  of  the  circular  hole  have  involved 
monotonic  uniaxial  remote  loading.  The  purpose  of  this  section  is  to 
study  the  effect  of  non-proportional  remote  loading  on  plastic  zone 
growth  and  compare  the  final  state  of  elastoplastic  deformation  for 
each  loading  path. 

The  finite  element  mesh  for  this  analysis  is  shown  in  Fig.  4.23. 

The  mesh  used  in  the  previous  analyses  was  modified  so  that  the 
external  dimensions  are  equal,  and  the  fixed  grip  end  condition 
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Figure  4.22  ANPLAST  residual  stress  distribution  after 
unloading  from  a  remote  stress  which  was 
twice  that  necessary  to  initiate  yielding. 
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Figure  4.22  continued 
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Figure  4.23  Finite  element  mesh  of  a  sheet  with  a  circular 
hole  for  analysis  of  complex  loading  paths. 


removed.  Three  separate  loading  cases  analyzed  are  as  follows 

I)  Vertical  Remote  Stress  to  41000  psi  followed  by  a 

Horizontal  Remote  Stress  to  41000  psi 

II)  Horizontal  Remote  Stress  (s^^  )  to  41000  psi  followed 
by  a  Vertical  Remote  Stress  i^yy)  to  41000  psi 

III)  A  simultaneous  application  of  both  remote  stresses 
^xx  ^yy  ^1000  Psi- 

The  isotropic  elastoplastic  material  properties  from  the  earlier 
analyses  are  used  in  all  three  cases. 

Discussion  of  Plastic  Zone  Growth 

The  plastic  zone  growth  for  Cases  I  and  II  illustrated  in 
Figs.  4.24-4.25  are  identical  except  that  they  are  oriented  oppositely 
about  the  respective  symmetry  axes.  In  both  of  these  cases,  the 
plastic  zone  progression  during  the  uniaxial  portion  of  the  loading  is 
of  the  same  character  that  was  observed  in  Fig. 4. 7  fur  the  fixed  grip 
analyses.  After  achieving  a  uniaxial  remote  load  of  41000  psi,  the 
opposing  transverse  remote  load  is  applied  incrementally  with  interest¬ 
ing  results.  All  but  six  elements  unload  into  the  elastic  regime, 
until  the  transverse  load  reaches  33000  psi  at  which  time  additional 
yielding  occurs  at  the  hole  perimeter  and  remote  locations. 

The  resulting  re-yielding  of  material  near  the  hole  does  not  take 
on  the  shape  one  would  expect  for  a  monotonically  increasing  uniaxial 
load.  The  stress  state  is  more  complex  due  to  the  biaxial  nature  of 
the  remote  loads,  and  the  previous  work-hardening  has  altered  the 
respective  effective  yield  stress  for  all  of  the  previously  yielded 
material.  Further  complications  result  from  the  assumed  non-linear 
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Figure  4.24  Plastic  zone  growth  for  an  isotropic  sheet  subjected  to  non-proportional 
remote  loading  (Case  I). 
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th  for  an  isotropic  sheet  subjected  to  non-proportional 


work-hardening  law.  Each  element  which  has  been  plastically  deformed 
can  be  at  a  different  location  on  the  effective  stress-effective 
plastic  strain  relation. 

The  plastic  zone  visualization  for  Case  III  shown  in  Fig. 4. 26 
represents  a  simultaneous  application  of  the  remote  loads.  Although 
the  loading  is  biaxial,  the  remote  components  are  stressed  at 
identical  rates  which  results  in  proportional  loading.  Both  the 
geometry  and  loading  are  symmetric  at  all  times  as  is  the  plastic  zone 
progression.  With  most  of  the  shear  stresses  removed  due  to  the 
equibiaxial  loading,  the  plastic  zone  progresses  in  rings.  It  is 
interesting  to  note  that  unlike  the  previous  two  cases,  plastic 
deformation  never  reaches  any  of  the  remote  loading  boundaries. 

Comparison  of  Permanent  Deformation  at  the  Final  Loading  Level 

The  previous  study  of  plastic  zone  f.  egression  clearly  demonstra¬ 
tes  that  the  distribution  cT  currently  yielded  material  is  different 
for  each  loading  path  r’ven  though  the  final  remote  stress  state  is 
the  same.  This  finding  suggests  that  plastic  deformation  in  each 
case  is  likely  to  be  different.  A  shortcoming  of  this  particular 
plastic  zone  visualization  technique  is  its  inability  to  relay 
quantitative  information  regarding  the  level  of  plastic  straining. 
Whether  an  element  has  experienced  extensive  plastic  deformation  or 
has  just  yielded,  the  same  asterisk  (*)  is  placed  at  the  centroid. 

An  iso-value  contour  of  the  plastic  deformation  would  provide 
the  necessary  quantitative  information.  A  suitable  measure  of  plastic 
deformation  for  this  purpose  is  the  effective  plastic  strain  (e^). 
Contours  of  e^  for  Cases  I,  II  and  III  are  presented  in  Figs.  4.27. 
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Figure  4.26  Plastic  zone  growth  for  an  isotropic  sheet  subjected 
to  a  simultaneous  biaxial  remote  loading  (Case  III). 
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It  is  obvious  from  these  contours  that  the  final  distribution  and 
magnitudes  of  permanent  deformation  anj quite  different  for  each  loading 
path,  even  though  the  final  remote  stress  states  are  identical.  These 
findings  substantiate  the  premise  of  load  path  dependency  and  highlight 
the  necessity  of  an  incremental  plasticity  approach. 

COMPARISON  TO  EXPERIMENTAL  RESULTS  FOR  A  CONTINUOUS  FIBER 
REINFORCED  METAL  MATRIX  COMPOSITE 

An  experimental  study  of  the  elastoplastic  response  of  a  unidirec¬ 
tional  Boron-Aluminum  metal  matrix  composite  strip  with  a  circular 
hole  was  conducted  by  S.  Rizzi  [24].  This  particular  experiment 
involved  loading  a  B-Al  specimen  as  shown  in  Fig.  A.ijs.  Multiple 
strain  gage  recordings  were  taken  in  the  vicinity  of  the  hole.  The 
finite  element  mesh  used  to  analyze  the  test  specimen  is  the  same  as 
for  the  other  fixed  grip  analyses  except  the  dimensions  were  scaled  to 
those  of  Fig.  4.28. 

The  test  was  run  by  increasing  the  load  until  certain  present 
breakpoints  were  reached  at  which  time  the  specimen  was  unloaded  to 
zero.  Since  most  of  the  strains  were  recorded  continuously,  a  strain 
distribution  at  any  load  can  be  easily  obtained. 

An  elastic  analysis  was  conducted  to  reconfirm  ANPLAST's  predictive 
capability  for  an  orthotropic  material.  The  elastic  properties  used 
are  those  described  for  B-Al  in  Reference  [12],  The  test  specimen  did 
not  appear  to  demonstrate  any  permanent  deformation  until  arriving  at 
a  remote  stress  of  approximately  4000  psi .  Therefore,  to  ensure  a 
purely  elastic  response  the  comparison  is  made  at  a  much  lower  load. 

The  experimental  results  shown  in  Fig.  4.29  are  plotted  against  the 
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Figure  4.28  Geometry  of  test  specimen  used  by  Rizzi  [24]  for 
experimental  determination  of  the  stress-strain 
behavior  of  a  Boron-Aluminum  metal  matrix  composite 
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Figure 


4.29  Comparison  of  ANPIAST  strains  with  the  experimental 

results  of  Rizzi  [24]  for  a  B-Al  sheet  with  a  circular 
hole  subjected  to  uniaxial  elastic  loading. 
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finite  element  solution  at  a  remote  stress  of  1700  psi.  The 
correlation  between  the  two  solutions  is  acceptable,  the  differences 
attributable  mainly  to  the  inherent  inaccuracies  of  measuring 
low  strain  levels. 

The  elastoplastic  material  properties  were  also  chosen  from 
Reference  [24].  The  elastic-plastic  loading  pathwas  chosen  so  that  it 
follows  the  test  procedure  as  close  as  possible.  The  elastic  analysis 
revealed  that  initial  yielding  would  occur  at  the  horizontal  hole 
perimeter  (e  *  0°)  at  a  remote  stress  of  4000  psi.  Therefore,  the 
first  loading  step,  which  must  be  elastic,  goes  straight  to  4000  psi, 
and  then  the  remainder  of  the  elastic  plastic  loading  steps  vary  in 
magnitude  but  are  never  larger  than  500  psi.  The  ANPLAST  loading  path 
explicitly  defined  in  Fig.  4.30.  Notice  that  increments  2,  6,  14  and 
17  coincide  with  the  test  unload  points.  The  desired  residual  stresses 
and  strains  are  obtained  by  re-starting  the  analyses  at  these  unload 
points  for  elastic  unloading.  Re-yielding  of  elements  during  unloading 
was  never  observed. 

The  plastic  zone  visualization  in  Fig.  4.31  follows  the  expected 
pattern  developed  throughout  this  section.  At  the  maxiimim  remote 
stress  of  10700  psi,  the  plastic  zone  has  progressed  completely  across 
the  specimen.  Along  the  horizontal  symmetry  line,  plasticity  has  also 
spread  across  the  net  section,  but  only  during  the  last  couple  loading 
increments.  Therefore,  it  is  reasonable  to  assume  that  little  or  no 
permanent  strain  will  be  found  at  the  outer  strain  gage  locations.  To 
evaluate  the  analysis  more  critically,  a  direct  study  of  the  elastic- 
plastic  strains  is  in  order. 


Figure  4.30  ANPLAST  Incremental  loading  path  for  analysis 
of  the  B-Al  strip  with  a  circular  hole. 


Figure  4.31  Plastic  zone  growth  for  the  B-Al  strip  subjected 
to  elastic-plastic  loading. 


The  longitudinal  and  lateral  strain  distributions  for  three 
elastic-plastic  remote  stresses  (6000  psi,  7700  psi,  10700  psi)  are 
plotted  against  the  experimental  results  in  Figures  4.32- 
4.33.  The  residual  strains  produced  after  unloading  are  included 
within  each  plot.  The  agreement  between  the  numerical  solution  and 
the  experimental  results  are  quite  good. 

Due  to  the  very  stiff  nature  of  the  reinforcing  fibers,  the 
lateral  strains  are  small  in  magnitude  and  both  the  analysis 

and  experiment  demonstrate  this.  The  longitudinal  (Cyy)  residual 
strains  demonstrate  that  ANPLAST  does  a  good  job  of  predicting  the 
plastic  deformation.  The  lack  of  any  substantial  residual  strains  at 
the  remote  gage  locations  is  in  agreement  with  the  plastic  zone 
visualization.  It  is  interesting  to  note  that  even  though  net  section 
yielding  occurs,  the  elastic  strains  still  dominate,  even  at  the  hole. 

Figure  4.34  highlights  the  stress-strain  response  at  the  innermost 
strain  gage  during  the  complete  load-time  history  of  the  test.  These 
results  substantiate  the  premise  that  ANPLAST  more  than  adequately 
predicts  the  response  of  structures  which  exhibit  anisotropic 
elastoplastic  behavior. 
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Figure  4.32  Comparison  of  ANPLAST  strains  with  experimental  results 
of  Rizzi  [24]  for  the  B-Al  strip  with  a  circular  hole 
subjected  to  elastic-plastic  loading. 
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Figure  4.33  Comparison  of  ANPLAST  lateral  strains  with  the 
experimental  resuUs  of  Rizzi  [24]  presented 
on  expanded  scale. 
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Figure  4.34  Comparison  of  ANPLAST  longitudinal  strain  with  the 
experimental  results  of  Rizzi  [24]  for  the  complete 
load-history  of  the  test.  The  strain  gage  is 
located  at  the  innermost  position  along  the 
horizontal  symmetry  line. 
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SECTION  5  ELASTOPLASTIC  STUDY  OF  CRACKED  PANELS 

INTRODUCTION 

The  study  of  fracture  and  fatigue  crack  propagation  often  requires 
an  understanding  of  the  deformation  state  in  the  vicinity  of  the 
crack  tip.  The  purpose  of  this  section  is  to  investigate  some  of  the 
material  and  loading  parameters  which  effect  crack  tip  plasticity. 
Special  attention  will  be  given  to  the  relative  effect  of  anisotropy 
on  plastic  zone  growth. 

This  section  will  first  deal  with  the  elastic  study  of  center- 
cracked  and  edge-cracked  panels  using  ANPLAST.  Two  numerical 
techniques  for  extracting  the  Mode  I  stress  intensity  factor,  will  be 
presented.  An  elastoplastic  analysis  of  large  scale  yielding  in 
isotropic  and  orthctropic  edge  cracked  panels  will  also  be  conducted 
using  ANPLAST. 

ELASTIC  ANALYSIS  OF  CENTER  CRACKED  PANELS 
The  center  cracked  panel  provides  a  benchmark  problem  to  further 
verify  the  performance  of  ANPLAST.  This  geometry  also  furnishes  the 
opportunity  to  study  the  affect  of  anisotropy  on  the  stress  distribu¬ 
tions. 

Comparison  of  Stress  Distribution  with  Theory 

The  center  cracked  panel  provides  the  simplest  geometry  for 
consideration  of  the  elastic  stress  distribution  in  the  vicinity  of  the 
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Figure  5.1  Finite  element  mesh  of  a  center  cracked  panel 
with  a  uniaxial  remote  load  and  fixed  grip  end 
conditions. 
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Principal  Material 
T  Direction 


0^  =  11  Material  Direction 
0  =  Angular  Orientation  of  Element 
r  =  Radius  of  Element 


Figure  5.2  Illustration  of  the  local  and  global  coordinate 
systems  for  a  cracked  panel . 
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Figure  5.3  Comparison  of  the  ANPLAST  stress  distribution 
with  the  ELLPLAS  theoretical  solution  for  an 
isotropic  panel  with  a  center  crack  subjected 
to  elastic  loading. 
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Figure  S4  Comparison  of  the  ANPLAST  stress  distribution 
the  ELLPLAS  theoretical  solution  for  a  weakly 
orthotropic  (B-Al)  panel  with  a  center  crack 
subjected  to  elastic  loading. 
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Figure  5.5  Comparison  of  the  ANPLAST  stress  distribution  with 
the  ELLPLAS  theoretical  solution  for  a  strongly 
orthotropic  (G-Ep)  panel  with  a  center  crack 
subjected  to  elastic  loading. 
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crack.  A  finite  panel  with  a  center  crack  as  illustrated  in  Fig.  5.1, 
wasanalyzed  with  ANPLAST  and  the  results  compared  to  the  theoretical 
solution  provided  by  ELLPLAS.  Isotropic  and  anisotropic  materials  are 
investigated. 

In  the  following  discussions  notational  definitions  are  as  shown 
in  Fig.  5.2.  The  stress  distributions  along  a  horizontal  radial  line 
(e  =  0°)  are  presented  in  Figs.  5. 3-5. 5.  The  ANPLAST  distribution  of 
stresses  for  the  isotropic  panel  and  both  B-Al  orientations  are  in 
close  agreement  with  the  theoretical  solution.  The  more  strongly 
aniostropic  graphite-epoxy  panel  exhibits  different  behavior, 
especially  when  the  fibers  are  parallel  to  the  crack  (e^  =  0°).  In 
this  case,  dominates  In  the  general  vicinity  of  the  crack  tip. 

The  numerical  model  of  the  center  cracked  panel  represents  a 
finite  width  and  length  sheet.  The  ELLPLAS  solution  is  that  for  an 
infinite  sheet.  Though  the  finite  effect  is  minimal  for  the  isotropic 
material  and  weakly  anisotropic  B-Al,  it  appears  that  strong  anisotropy 
can  amplify  the  finiteness  of  the  specimen.  This  is  evident  from  the 
behavior  of  in  the  graphite  epoxy  center  cracked  sheet.  When  the 
fibers  are  parallel  to  the  crack,  decays  dramatically  over  the 
width  of  the  specimen,  while  for  =  90°,  this  stress  component  has 
dropped  to  almost  zero  a  very  small  distance  from  the  crack  tip.  These 
results  indicate  that  if  an  accurate  elastic  analysis  is  desired  in 
the  vicinity  of  the  crack  tip,  mesh  density  should  be  a  factor  of 
anisotropy  as  well  as  geometry. 


125 


Philosophy  of  Finite  Element  Mesh  Design 

The  analysis  of  the  center  cracked  panel  from  the  previous  section 
used  a  radially  designed  finite  element  mesh.  A  symmetric  radial  fan 
is  centered  about  the  crack  tip  and  the  mesh  density  is  biased  toward 
the  crack  tip.  An  alternative  to  this  mesh  design  would  be  a  rectang¬ 
ular  mesh  as  illustrated  in  Fig.  5.6.  One  possible  advantage  for 
selecting  the  rectangular  design  would  be  the  ease  of  mesh  generation. 

A  concern,  is  the  relative  accuracy  in  representing  the  stress 
distribution  in  the  vicinity  of  the  crack. 

An  elastic  analysis  was  conducted  utilizing  the  rectangular  mesh; 
the  results  are  presented  in  Fig.  5.7.  The  numerical  results  are  again 
plotted  against  the  theoretical  solution  and  compared  to  the  respective 
plot  for  the  radially  designed  mesh.  While  the  stress  distributions 
for  the  rectangular  mesh  are  reasonable,  it  is  apparent  for  e  *  90°, 
that  the  radial  design  is  superior.  It  is  interesting  to  note  that  the 
rectangular  mesh  contains  104  more  d.o.f.  than  does  the  radial  mesh. 

An  elastic  analysis  was  used  to  compare  the  two  mesh  designs,  yet 
it  is  the  study  of  plasticity  which  is  given  primary  emphasis  in  this 
research.  When  considering  the  Initiation  of  yield  and  subsequent 
plastic  flow,  one  must  recall  that  these  processes  are  calculated  on 
the  element  level.  It  is  obvious  that  the  radial  mesh  design  provides 
many  more  elements  in  the  e  direction,  thus  allowing  a  better 
description  of  the  gradients  associated  with  plasticity.  Therefore, 
all  elastoplastic  analyses  of  cracked  panels  will  utilize  a  radially 
designed  mesh. 
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Figure  5.6  Finite  element  mesh  of  a  center  cracked  panel  which 
utilizes  a  rectangular  grid  design  in  the  vicinity 
of  the  crack  tip. 
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Figure  5.7  Comparison  of  ANPLAST  stress  distribution  for  the 

radial  and  rectangular  mesh  designs  of  an  isotropic 
center  cracked  panel  subjected  to  elastic  loading. 
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Comparison  of  Constant  Strain  Triangle  and  Isoparametric  Quadrilateral 
Elements 

The  primary  finite  element  used  in  the  ANPLAST  program  is  the 
constant  strain  triangle  (CST).  It  has  been  suggested  by  a  number  of 
authors  [23,  25]  that  the  higher  order  isoparametric  element  with  a 
reduced  integration  scheme  may  provide  superior  results  for  analyses 
involving  strong  anisotropy.  Therefore,  the  8  node  isoparametric 
element  developed  in  Section  3  was  used  to  analyze  the  center  cracked 
panel  and  the  results  compared  to  the  CST.  The  mesh  used  for  this 
analysis  is  illustrated  in  Fig.  5.8. 

The  results  from  the  comparative  study  are  pi  esented  in  Fig.  5.9. 
The  results  are  presented  for  an  isotropic,  B-Al  (e^  ®  0°  ),  and  6-Ep 
(e^  =  90°)  material.  Note  the  erratic  behavior  o1  the  two  nodal 
stresses  closest  to  the  crack  tip.  This  is  due  tc  the  nodal  extra¬ 
polation  technique  failing  to  account  for  the  8  node  quadralateral 
being  degenerated  to  a  6-node  triangle. 

Of  interest  is  the  0^^^^  distribution  for  the  graphite  epoxy 
specimen  since  this  resulted  in  the  greatest  deviation  for  the  CST.  It 
appears  that  the  quadrilateral  element  does  provide  a  better  represen¬ 
tation  of  this  stress  gradient.  But  rather  than  the  8  noded  quadri¬ 
lateral  possessing  any  inherent  advantage,  the  improved  performance 
is  likely  due  to  the  available  degrees  of  freedom.  The  quadrilateral 
mesh  has  almost  twice  the  degrees  of  freedom  than  does  the  CST  mesh. 
Therefore,  a  more  refined  CST  mesh  should  produce  results  similar  to 
the  8  node  quadrilateral,  for  approximately  the  same  cost  (running 
time).  The  use  of  either  element  in  a  appropriately  refined  mesh 
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Figure  5.9  Comparison  of  the  Isoparametric  element  results  with 
the  ELLPLAS  solution  for  the  center  cracked  panel 
subjected  to  elastic  loading. 
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Figure  5.9  continued 
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should  result  in  adequate  performance. 


ELASTIC  ANALYSIS  OF  AN  EDGE  CRACKED  PANEL 

All  of  the  material  behavior  investigated  in  this  research  has 
exhibited  an  elastic  response  before  the  onset  of  plasticity,  and  it 
is  clear  that  the  elastic  stress  distribution  plays  a  significant 
role  in  subsequent  plasticity.  In  addition,  the  elastic  behavior  of 
cracks  is  analyzed  because  of  the  necessity  for  understanding  some 
basic  fracture  mechanics  parameters  and  how  they  are  influenced  by 
anisotropy. 

Effect  of  Anisotropy  on  Stress  Distribution 

It  is  reasonable  to  expect  that  the  elastic  stress  distribution 
in  the  vicinity  of  the  crack  tip  for  the  edge  cracked  panel  retains 
the  same  character  as  was  discussed  for  the  center  crack.  For  cases 
of  isotropic  or  weakly  anisotropic  material  behavior,  will  dominate 
along  the  horizontal  cross  section.  The  introduction  of  strong 
anisotropy  will  exaggerate  the  finiteness  of  the  panel  and  redistribute 
the  stress  gradients.  Therefore,  similar  modeling  techniques  were 
utilized  for  analyses  of  the  edge  cracked  panel  as  was  tsed  for  the 
center  crack. 

A  detailed  understanding  of  the  individual  stress  components  does 
not  provide  any  direct  information  concerning  the  onset  of  yield.  But 
an  equivalent  quantity  such  as  the  Von>Mises  [2]  effective  stress 
parameter  provides  useful  insight  into  isotropic  elastoplasticity.  A 
parametric  study  on  the  effect  of  elastic  anisotropy  via  the  Von-Mises 
stress  can  be  conducted  in  a  meaningful  way  only  if  one  assumes  all 
plasticity  behavior  is  isotropic. 


Five  elastic  analyses  of  an  edge  cracked  panel  as  shown  in  Fig. 
5.10  were  conducted  for  an  isotropic,  B-Al  (e^  =  0°,  90°),  and  G-Ep 
(e^  =  0°,  90°)  material.  The  results  are  in  the  form  of  iso-stress 
contours  of  the  Von-Mises  quantity  as  illustrated  in  Figs.  5.11-5.13. 
Since  all  of  the  materials  were  assumed  to  behave  in  an  isotropic 
plastic  manner,  the  Von-Mises  contours  leads  us  to  the  conclusion  that 
plastic  zone  shape  is  strongly  linked  to  the  strength  and  orientation 
of  the  elastic  anisotropy.  This  is  clearly  illustrated  for  the  G-Ep 
specimen  response,  which  demonstrates  the  significant  distortion  of 
the  contours  in  the  stiff  material  direction.  By  noting  the  maximum 
contour  level,  one  can  conclude  that  the  remote  load  which  initiates 
yield  is  also  a  function  of  the  strength  and  orientation  of  the 
elastic  anisotropy. 

It  is  clear  that  elastic  anisotropy  can  have  a  substantial  effect 
on  the  elastoplastic  response.  Earlier  sections  have  documented  the 
consequences  of  introducing  anisotropy  into  the  plastic  flow  relation¬ 
ships.  For  orthotropic-elastic  orthotropic- plastic  materials,  these 
two  effects  will  superimpose  and  result  in  an  extremely  complex 
structural  response.  This  highlights  the  convenience  of  using  a 
numerical  tool,  like  ANPLAST,  to  analyze  these  types  of  problems. 

Stress  Intensity  Determination  by  an  Energy  Technique 

The  advantage  of  using  an  energy  technique  to  calculate  fracture 
parameters  lies  in  the  fact  that  no  extrapolation  or  curve-fitting  is 
needed.  In  particular,  the  strain  energy  release  rate,  g^,  of  a 
cracked  panel  can  be  derived  by  investigating  the  change  in  compliance 
of  the  system  as  the  crack  extends.  It  can  be  shown  [26]  that  the 
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Figure  5.10  Finite  element  mesh  of  a  transverse  edge  cracked 
panel  with  a  uniaxial  remote  load  and  fixed  grip 
end  conditions. 
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Figure  5.11  Mises  stress  contour  of  a  isotropic  edge  cracked 
panel  subjected  to  elastic  loading. 
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panel  subjected  to  elastic  loading. 
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panel  subjected  to  elastic  loading. 


change  in  elastic  strain  energy  during  crack  extension  is  equal  to  the 
change  in  potential  energy  (aV)  and  can  be  expressed  as 

aV  »  -[  [  Pd6  -  f  Pd6]  *  -  i  PCup-Uj)  (5.1) 

Jq  Jq 

where  (U2-Uj)  is  the  deflection  of  the  system  in  the  direction  (6)  of 
the  applied  force  (P),  at  the  point  where  the  force  is  applied.  The 
change  in  strain  energy  (aU)  can  now  be  equated  to  the  compliance  of 
the  system 


AU  «  AV  «  ^  P(Au)  (5.2) 

and  the  strain  energy  release  rate  defined  as  the  change  in  strain 
energy  due  to  a  small  increment  of  crack  extension  (Aa)  becomes 
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(5.3) 


It  is  now  clear  that  gj  can  be  estimated  by  measuring  the  deflections 
of  the  system  at  the  point  of  load  application. 

The  compliance  method  is  ideally  suited  for  use  with  the  finite 
element  method.  Since  the  finite  element  technique  directly  solves  for 
displacements,  no  interpolation  of  stresses  or  strains  is  required. 

The  work  done  by  external  forces  is  easily  calculated  by  multiplying 
the  applied  nodal  loads  by  the  respective  change  in  displacements  due 
to  a  small  crack  extension.  It  is  interesting  to  note  that  deriving  gj 
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by  the  compliance  method  does  not  require  any  special  technique  or 
program  when  investigating  orthotropic  materials. 


Although  g^  can  be  used  directly  in  many  fracture  and  fatigue 
crack  propagation  applications,  it  is  still  quite  common  to  utilize 
the  stress  intensity  factor  (Kj).  The  relation  between  gj  and  Kj  for 
an  isotropic  material  under  a  plane  stress  condition  is  [26] 


9l 


(5.4) 


where  E  is  the  modulus  of  elasticity.  After  calculating  gj,  it  is 
easy  to  find  Kj.  For  an  orthotropic  material,  the  relation  is  not 
readily  available,  but  is  in  the  form 


where 


6  *  f(Ej,  ^2*  '*12*  ^12^  (^*5) 

To  determine  6  analytically  is  cumbersome,  but  taking  advantage  of  the 
properties  of  the  center-cracked  panel,  B  can  be  found  numerically. 

For  a  center-cracked  panel,  where  the  finite  width  effects  are 
negligible,  Kj  is  identical  for  both  the  isotropic  and  orthotropic 
material  cases  [27],  Therefore,  Kj  can  be  found  from  the  simple 
expression 
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Using  the  compliance  method  previously  explained,  and  the  center- 
cracked  finite  element  model,  gj  is  calculated  for  a  B-Al  specimen 
with  the  fibers  oriented  parallel  to  the  crack.  Since  B  was  calculated 
and  found  to  have  a  value  of  9.11x10^  psi,  it  is  possible  to  find 
for  any  B-Al  edge-cracked  geometry. 

A  stress  intensity  calibration  chart  prepared  for  an  edge-cracked 
panel  is  illustrated  in  Fig.  5.14.  The  solid  line  represents  the 
solution  for  an  isotropic  material,  with  bending  restrained.  The  four 
finite  element  meshes  used  for  these  analyses  are  similar  to  the  mesh 
in  Fig.  5.10  except  for  the  length  of  the  crack.  The  triangular 
data  points  represent  Kj  derived  by  extrapolation  of  the  stresses 
as  discussed  in  the  next  section. 

Stress  Intensity  Determination  by  Extrapolation  of  the  Stresses 

The  determination  of  the  Mode  I  stress  intensity  factor  Kj  by 
direct  extrapolation  of  the  stresses  is  of  interest  in  this  section. 

For  anisotropic  materials,  the  near  field  stresses  can  be  represented 
in  terms  of  Kj,  r,  and  9  as  [27] 


,  ^  f^^(9) 
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(5.6) 
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Figure  5.14  Stress  Intensity  Factor  (Kj)  chart  for  the  edge 

cracked  panel  containing  the  isotropic  handbook 
solution  and  that  obtained  by  ANPLAST  for  B-A1. 
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where  *  appropriate  functional  form  of  e.  This  equation  can 

be  easily  re-arranged  so  that  Kj  appears  on  the  left  hands ide  of  each 
expression  as 


o  .  .  /r  /2it 

(5-7) 

Assuming  that  is  known,  it  is  apparent  that  Kj  can  be  obtained 

from  a  finite  element  analysis  since  the  stresses  o.^  at  each  location 
(r,  0)  are  known.  Ideally,  any  selection  of  stress  points  would 
result  in  the  same  value  of  Kj.  Unfortunately,  if  stresses  are 
sampled  from  the  FEM  results  at  various  angular  and  radial  locations, 
differing  values  of  Kj  are  obtained.  Because  of  this  variation,  an 
extrapolation  technique  along  with  a  least  squares  curve  fitting 
procedure  must  be  used. 

The  first  step  is  to  design  a  finite  element  mesh  which  centers 

a  radial  fan  on  the  crack  tip.  This  mesh  orientation  provides  many 

lines  of  constant  0  on  which  the  extrapolation  technique  can  be  used. 

After  the  finite  element  analysis  is  completed,  the  nodal  stresses 

are  post-processed  and  plotted  along  radial  lines,  as  shown  in  Fig. 

5.15.  Note  that  in  this  figure  the  ordinate  is  (a  +  o  )  /r  and 

XX  yy 

the  abscissa  is  simply  /?  .  It  was  found  in  Reference  [16]  that  the 
sum  of  the  normal  stress  components  behave  in  a  more  predictable 
manner,  so  this  technique  was  utilized  here.  The  reason  for  plotting 
the  stresses  in  this  manner  can  be  demonstrated  by  again  rearranging 
Eq.  5.6  written  as 
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Figure  5.15  Stress  distribution  along  19  radial  lines 

emanating  from  the  crack  tip  of  the  transverse 
edge  cracked  panel. 
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Since  Kj  is  defined  in  the  limit  as  the  radial  distance  approaches 
zero,  one  can  extrapolate  the  results  in  Fig.  5.15  for  each  radial 
line.  In  essence,  each  line  is  being  extended  to  finds  its  intercept  with 
the  y-axis.  This  is  done  in  a  least  squares  sense  by  curve  fitting 
a  straight  line  through  each  set  of  data  points.  Usually,  the  nodal 
stress  closest  to  the  crack-tip  is  excluded  since  it  appears  to  be 
distorted.  Now,  for  each  of  the  nineteen  (19)  radial  orientations, 

(o  +  0  )  /r  is  obtained. 

AA  yjr 

The  next  step  is  to  use  ELLPLAS  to  provide  f^^(e).  The  advantage 
of  utilizing  ELLPLAS  is  that  it  can  furnish  fi^(e)  for  art  orthotropic 

’  w 

material  as  well  as  an  isotropic  material.  Entering  some  value  of 
and  multiplying  by  f^.(e),  the  left  hand  side  of  Eq.  5.8  can  be 
calculated.  is  now  determined  by  using  a  least  squares  fit  of  the 
data  provided  in  the  previous  exercise.  The  results  for  an  isotropic 
and  B-Al  edge-crack  are  illustrated  in  Fig.  5.16.  It  is  interesting 
to  note  that  the  stress  intensity  factor  Kj  extrapolated  for  both 
cases  was  the  same  (17000  psi  /in)  and  is  probably  due  to  the  weak 
elastic  anisotropy  exhibited  by  B-Al. 

It  is  apparent  from  the  results  presented  that  the  energy  method 
is  the  most  accurate  and  simplest  for  finding  Kj.  Although  it  requires 
two  analyses  to  be  run,  the  post-processing  of  data  is  minimal  and 
the  results  are  very  good.  The  stress  component  extrapolation  technique 
was  investigated  as  a  means  of  determining  the  quality  of  the  field 
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Figure  5.16  Stress  Intensity  Factor  calibration  plot.  A 
least  squares  curve  fit  of  Kj  F^j(e)  vs  the 

19  radial  intercept  values 
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quantities  for  subsequent  plasticity  analysis.  In  general,  the 
quality  of  these  results  is  good. 

ISOTROPIC  ELASTIC-PLASTIC  ANALYSIS  OF  AN  EDGE  CRACK 

The  isotropic  elastoplastic  behavior  of  an  edge  cracked  panel  was 
investigated  to  provide  a  benchmark  solution  from  which  orthotropic 
elastoplastic  response  can  be  compared. 

A  Study  of  Plastic  Zone  Growth 

An  isotropic  elastic-plastic  analysis  of  the  transverse  edge- 
crack  was  conducted.  Work-hardening  was  expressed  by  the  Ramberg-Osgood 
relationship  discussed  in  the  previous  sections.  The  progression  of 
the  plastic  zone  is  illustrated  in  Fig.  5.17.  Note  the  distinctive 
lobe  at  about  60°,  where  the  Von-Mises  yield  criterion  [2  ]  predicts 
the  maximum  extent  of  the  plastic  zone.  Fig.  5.18  provides  additional 
insight  into  the  effective  stress  distribution  and  how  it  ultimately 
controls  the  plastic  zone  development.  At  a  remote  stress  of  24,500 
psi,  the  plastic  zone  has  progressed  more  than  halfway  across  the 
width  of  the  specimen,  indicating  that  large  scale  yielding  is  being 
considered. 

The  consequences  of  mixed  mode  loading  on  plastic  zone  growth 
investigated  by  the  analysis  of  an  oblique  crack  illustrated  in  Fig. 

5.19.  The  remote  loading  is  uniaxial,  but  shear  stresses  are  now 
introduced  along  the  plane  of  the  crack.  The  plastic  zone  progression 
illustrated  in  Fig.  5.20  demonstrates  a  significant  departure  from  the 
transverse  crack  behavior.  The  effective  stress  contours  in  Fig.  5.21 

substantiate  the  qualitative  plastic  zone  representation.  Considering  the 
previously  discussed  analogy  which  relates  plastically  deformed 
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Figure  5.17  Plastic  zone  growth  for  an  isotropic  edge-cracked 
panel  subjected  to  elastic-plastic  loading. 
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Figure  5.18  Effective  stress  contour  of  an  isotropic  edge^cracked 
panel  at  the  onset  of  yielding. 
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Figure  5.19  A  finite  element  mesh  of  an  oblique  edge-cracked 
panel  subjected  to  uniaxial  loading  and  fixed 
grip  end  conditions. 
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Figure  5.20  Plastic  zone  growth  for  an  isotropic  oblique 
edge-cracked  panel  subjected  to  elastic- 
plastic  loading. 
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Figure  5.21  Effective  stress  contours  for  the  oblique  edge-cracked  panel 
at  the  onset  of  yielding. 


material  to  damage  initiation,  it  is  clear  that  under  cyclic  loading 
the  oblique  crack  will  tend  to  turn  and  propagate  horizontally,  which 
has  been  observed  in  the  laboratory  [28]. 

The  Effect  of  a  Remotely  Applied  Transverse  Stress 

In  most  practical  applications  the  remote  loading  is  rarely 
uniaxial.  For  an  elastic  analysis  of  a  crack  perpendicular  to  a 
remote  uniaxial  load,  classical  LEFM  suggests  that  the  addition  of  a 
remote  load  parallel  to  the  crack  will  not  contribute  to  Kj.  In  a 
superficial  manner,  this  suggests  that  a  transverse  biaxial  loading 
will  have  no  effect  on  fracture  or  fatigue  crack  propagation.  If 
plasticity  is  considered,  then  a  transverse  remote  stress  will  change 
the  character  of  the  problem. 

Three  different  ratios  (x)  of  transverse  to  axial  load  were  applied 
to  the  transverse  edge  cracked  model.  The  axial  l^yy)  remote  stress 
•«s  used  to  gauge  the  load  incrementation  and  the  transverse  stress 
(s^^)  was  added  simultaneously  in  the  stated  proportion.  Three  sets 
of  plastic  zones  progressions  are  illustrated  in  Fig.  5.22. 

The  lowest  level  of  transverse  loading  (x  *  0.5)  results  in  a 
change  in  both  size  and  shape  of  the  plastic  zone  in  comparison  to  the 
ourely  uniaxial  remote  loading.  The  effect  of  the  remote  biaxial 
stress  appears  to  be  a  flattening  of  the  60°  lobe.  At  a  remote  load 
of  7500  psi,  the  plastic  zone  for  the  uniaxial  case  is  more  than 
halfway  across  the  width  of  the  net  cross-section,  but  for  the 
transverse  loading  cases  (x  *0.5  and  x  *  1-0)  it  is  visibly  reduced. 
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Figure  5.22  Plastic  zone  progression  for  an  isotropic 
transverse  edge-cracked  panel  subjected  to 
proportional  biaxial  remote  loading. 
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Figure  5.22  continued 
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This  behavior  is  related  to  the  shape  of  the  isotropic  yield 
locus.  If  yield  has  been  reached  for  a  purely  uniaxial  stress,  the 
application  of  a  transverse  load  results  in  a  translation  inside  the 
yield  surface,  and  an  unloading  to  an  elastic  state.  For  the  edge 
crack  analysis,  the  biaxial  loading  is  proportional,  but  a  similar 
reduction  in  yielding  is  observed  with  respect  to  the  uniaxial  remote 
load. 

The  maximum  transverse  remote  loading  (x  =  2.0)  demonstrates 
significantly  different  plastic  zone  progression.  The  transverse  load 
appears  to  have  completely  dominated  the  Mode  I  component  of  the 
plastic  zone.  The  crack  tip  acts  as  a  concentration  which  initiates 
yielding,  but  the  transverse  loading  quickly  induces  net  section 
yielding. 

It  has  been  demonstrated  that  the  plastic  zone  is  significantly 
altered  due  to  a  transverse  remote  load.  Therefore,  it  is  reasonable 
to  assume  that  transverse  loading  will  in  turn  have  an  impact  on 
fracture  and  fatigue  behavior. 

A  Study  of  Unloading  and  the  Resulting  Residual  Stresses 

The  transverse  edge  cracked  specimen  was  loaded  into  the  elastic- 
plastic  regime  to  a  remote  load  of  2ip00  psi  and  then  unloaded.  This 
unloading  was  conducted  in  an  incremental  manner  to  allow  the  yielded 
material  to  become  elastic  again.  Indeed,  after  a  very  small 
unloading  increment,  the  entire  panel  (all  elements)  becomes  elastic, 
until  the  last  unloading  increment  where  a  few  elements  in  the  vicinity 
of  the  crack  tip  re-yield. 
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The  elastic  material  which  surrounds  the  plastic  zone  compresses 
the  permanently  deformed  material  upon  unloading.  The  outcome  of  this 
compressing  action  results  in  residual  stresses,  especially  near  the 
crack  tip  where  the  plastic  strains  are  the  greatest.  Fig.  5.23 
presents  the  residual  stress  distribution  for  the  above  unloading  case 
and  a  1.35  overload  followed  by  unloading.  In  both  cases,  the 
classical  compressive  stress  is  present  in  the  near  vicinity  of 

the  crack  tip.  It  is  the  compressive  residual  stress  component  which 
is  thought  of  as  retarding  any  subsequent  fatigue  crack  propagation. 

ORTHOTROPIC  ELASTIC-PLASTIC  ANALYSIS  OF  AN  EDGE  CRACK 
The  effect  of  anisotropy  was  considered  in  the  analysis  of  the 
edge  cracked  panel  to  gain  a  better  understanding  of  the  differences 
from  isotropic  behavior. 

The  Effect  of  Anisotropy  on  Plastic  Zone  Growth 

The  progression  of  the  plastic  zone  in  an  orthotropic  (B-Al) 
elastoplastic  transverse  edge  cracked  panel  was  investigated.  The 
plastic  zone  visualization  is  presented  in  Fig.  5.24  for  the  case  where 
the  fibers  are  parallel  and  perpendicular  to  the  crack.  The  elastic 
effective  stress  contours  are  plotted  in  Fig.  5.25.  As  was  the  case 
for  the  isotropic  edge  crack,  large  scale  yielding  was  considered. 

The  results  for  the  transverse  crack  parallels  the  elastoplastic 
analysis  of  a  sheet  with  a  circular  hole.  The  plastic  zone  shape  for 
the  isotropic  material  and  B-AI  (e^  *  0°)  are  similar,  while  for 

*  90°  a  significant  departure  is  observed.  With  the  fibers  perpendicular 
to  the  crack  the  majority  of  plastic  deformation  occurs  in  a  vertical 
band  directly  above  the  crack  tip.  Unlike  the  circular  hole,  a  small 
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Figure  5.23  ANPLAST  residual  stress  distribution  following 


unloading  for  an  isotropic  transverse  edge 
cracked  panel. 
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Figure  5.24  Plastic  zone  growth  for  an  orthotropic  B-Al  transverse 
edge-cracked  panel  subjected  to  elastic-platic 
loading. 
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Figure  5.24  continued 
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Figure  5.25  Effective  stress  contours  of  an  orthotropic  B-Al  transverse 
edge-cracked  panel  at  the  onset  of  yielding. 


area  of  permanent  deformation  is  present  in  front  of  the  crack  tip. 

The  oblique  crack  mode!  was  also  analyzed  to  investigate 
plasticity  effects  in  B-Al  under  a  mixed  mode  loading.  Fig.  5 .26 
illustrates  the  growth  of  the  plastic  zone  for  the  case  where  the 
fibers  are  running  parallel  to  the  crack.  The  results  are  signifi¬ 
cantly  different  from  those  observed  for  the  isotropic  case.  The  B-Al 
plastic  zone  tends  to  produce  two  lobes:  one  directly  in  line  with 
the  crack  tip  (parallel  to  the  fibers)  and  another  at  120°  in  a 
clockwise  orientation. 

These  results  again  support  the  assumption  that  under  cyclic 
loading,  plastically  deformed  material  represents  a  likely  sight  for 
further  crack  initiation.  It  has  been  observed  in  the  laboratory  [12] 
for  B-Al  edge  cracked  specimens  that  fatigue  crack  propagation  proceeds 
parallel  to  the  fibers,  which  is  consistent  with  Fig.  5.26. 

The  Effect  of  a  Remotely  Applied  Transverse  Stress 

The  consequences  of  a  remotely  applied  transverse  stress  on  a 
B-Al  edge  cracked  panel  was  studied.  The  analysis  parallels  the  study 
for  the  isotropic  panel,  except  the  transverse  to  axial  load  ratios  (x) 
were  1.0,  2.0,  and  4.0  respectively.  The  plastic  zone  visualizations 
are  presented  in  Fig.  5.27.  Only  the  case  of  the  fibers  parallel  to 
the  crack  (e^  =  0°)  is  analyzed. 

The  plastic  zone  progression  for  all  load  ratios  demonstrates 
t^at  the  transverse  loading  has  almost  no  effect  on  plastic  deforma¬ 
tion.  Both  the  shape  and  relative  magnitude  of  these  plastic  zones 
are  unaltered.  This  behavior  can  be  explained  if  one  considers  the 
anisotropy  induced  by  the  fiber  orientation.  The  predominant  stress 
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Figure  5.26  Plastic  zone  growth  for  an  orthotropic  B-Al  oblique 
edge-cracked  panel  subjected  to  elastic-plastic 
loading. 
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Figure  5.27  Plastic  zone  growth  for  an  orthotropic  B-Al 
transverse  edge-cracked  panel  subjected  to 
proportional  biaxial  remote  loading. 
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Figure  5.27  continued 
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Induced  by  the  transverse  load,  is  multiplied  by  the 

plasticity  coefficient  rj^^  in  the  effective  stress  calculation.  Since 
r^^  is  only  0.001  for  B-Al,  the  transverse  stress  has  almost  no  effect 
on  yield. 

One  location  where  the  plastic  zone  does  appear  to  be  affected  is 
towards  the  fixed  grip  location.  Here,  the  transverse  stress  has  the 
effect  of  unloading  a  few  elements  due  to  the  tilt  of  the  yield  locus 
as  was  discussed  for  the  isotropic  material.  In  any  case,  one  can 
state  with  some  authority  that  transverse  loading  will  not  strongly 
influence  fracture  or  fatigue  of  a  B-Al  edge  cracked  specimen  with 
this  particular  material  orientation. 

A  Study  of  Unloading  and  the  Resulting  Residual  Stresses 

The  B-Al  edge  cracked  panel  was  subjected  to  two  separate  load 
excursions.  The  residual  stress  distributions  are  plotted  in  Fig. 

5.28.  The  results  exhibited  are  almost  identical  to  the  isotropic 
examples.  Rather  large  compressive  stresses  can  be  found  in  the 
vicinity  of  the  crack  tip,  with  dominating  in  both  cases. 

Again,  the  classical  line  of  thinking  is  that  residual  compressive 
stresses  tend  to  inhibit  further  fatigue  crack  propagation.  Therefore, 
one  can  conclude  that  any  deviation  from  classical  retardation  for 
orthotropic  plastic  materials  is  not  a  stress  related  phenomenon. 

A  Look  at  the  Crack  Opening  Profile 

It  has  been  suggested  [26]  that  the  crack  opening  displacement 
(COD)  can  be  used  as  a  nonlinear  fracture  mechanics  parameter  for 
elastic-plastic  materials.  It  is  assumed  that  the  displacement 
profile  of  the  crack  is  related  to  the  state  of  deformation  in  the 
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Figure  5.28  ANPLAST  residual  stress  distribution  following 
unloading  for  an  orthotropic  B-A1  transverse 
edge-cracked  panel. 
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vicinity  of  the  crack  tip,  yet  independent  of  the  particular  specimen 
geometry.  The  advantage  of  this  technique  is  obvious  to  the  experimen¬ 
talist.  The  crack  opening  displacement  is  relatively  easy  to 
measure  and  does  not  require  the  sophisticated  instrumentation 
necessary  when  trying  to  obtain  near  field  stresses  and  strains. 

The  crack  opening  profiles  for  the  isotropic  and  B-Al  {e^  =  0°, 
90°)  edge  cracked  analyses  are  presented  in  Fig. 5.29  .  Three  sets  of 
profiles  are  plotted;  elastic  loading,  elastic  plastic  loading,  and 
complete  unloading.  For  the  elastic  and  elastic-plastic  loading  cases, 
all  of  the  profiles  have  approximately  the  same  shape,  with  only  small 
deviations  occurring  a  small  distance  from  the  tip.  The  unloading 
profiles  demonstrate  a  much  stronger  contrast  in  shape.  The  peak 
value  of  displacement  is  different  with  respect  to  the  crack  opening. 
This  indicates  that  while  the  near  tip  displacement  is  approximately 
proportional  to  crack  opening  during  monotonic  loading,  the  amount  of 
elastic-plastic  near  tip  recovery  is  not. 

It  should  be  noted  that  the  triangular  fan  which  is  centered 
about  the  crack  tip  is  experiencing  a  large  amount  of  deformation  due 
to  the  discontinuity  in  geometry.  It  is  likely  that  a  single  linear 
displacement  element  is  not  adequate  to  represent  this  large  gradient. 
Therefore,  the  resulting  displacement  spike  is  not  realistic,  and  any 
errors  in  displacements  near  the  tip  '-*111  have  an  additive  effect. 

The  consequences  of  these  findings  have  not  been  fully  developed. 
It  may  turn  out  that  these  particular  profiles  can  be  related  to  a 
crack  closure  phenomenon  [29],  as  illustrated  by  the  contrast  of  the 
isotropic  and  B-Al  (e  •  0°)  profiles  after  unloading.  It  has  been 


171 


Figure 


5.29  Crack  opening  profiles  for  a  transverse  edge-cracked 
panel  under  uniaxial  remote  loading.  Note  that  the 
displacements  are  normalized  by  the  outer  edge 
displacement. 
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observed  in  the  laboratory  [12],  that  retardation  effects  in  B-Al 
(e  =  0°)  edge  cracked  panels  are  minimal  in  comparison  to  isotropic 
specimens  of  the  same  geometry.  It  could  be  proposed  that  the  lack 
of  a  "dipping"  of  the  B-Al  profile  negates  any  contact  with  the 
opposite  face  upon  unloading  which  could  produce  the  crack  closure 
effect.  But,  before  any  firm  conclusions  could  be  proposed,  additional 
research  would  be  required. 

Comparison  to  Experimental  Results  for  a  Continuous  Fiber  Reinforced 
Metal  Matrix  Composite 

An  experimental  study  of  the  elastoplastic  response  of  a 
unidirectional  Boron>Aluminum  metal  matrix  composite  edge  cracked  panel 
was  conducted  by  0.  Kenaga  [12]  and  S.  Rizzi  [24].  This  particular 
study  involved  testing  both  a  transverse  and  oblique  edge  cracked 
specimen,  as  illustrated  in  Fig.  5.30.  Multiple  strain  gage  recordings 
were  taken  in  the  vicinity  of  the  crack  tip.  The  finite  element 
meshes  used  to  analyze  these  test  specimens  are  the  same  as  for  the 
other  fixed  grip  analyses,  except  the  dimensions  were  scaled  to  those 
of  Fig.  5.30. 

The  tests  were  run  by  first  cycling  the  load  to  produce  a  distinct 
crack,  then  applying  the  strain  gages.  The  specimens  were  cycled 
again  until  the  crack  tip  propagated  to  the  vicinity  of  the  first 
strain  gage.  The  strain  gages  were  then  zeroed  and  a  monotonic 
overload  applied.  The  ANPLAST  analyses  incorporates  monotonic  loading 
to  the  maximum  stress. 

The  results  for  the  gage  nearest  the  crack  tip  along  the  plane 
of  the  crack  was  presented  in  Fig.  5.31.  In  both  cases  the  finite 
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Figure  5.31  Comparison  of  ANPLAST  strains  with  the  experimental  results  of  Kenaga  [12] 
and  Rizzi  [24]  for  a  B-Al  panel  with  a  transverse  and  oblique  edge  crack. 


element  results  overestimate  the  amount  of  work-hardening  measured.  It 
is  apparent  that  the  cyclic  loading  prior  to  the  overload  changed  the 
elastoplastic  material  properties  of  the  B-Al.  Obviously,  some  type 
of  damage  has  been  induced  into  the  B-Al  microstructure  by  the  cycling 
of  the  load  which  softens  the  resulting  plastic  work-hardening. 

The  material  properties  used  in  ANPLAST  to  analyze  these  problems 
were  derived  from  monotonically  loaded  uniaxial  specimens  with  no  prior 
load  conditioning.  If  a  more  representative  analysis  were  required, 
then  the  appropriate  elastic  and  plastic  material  parameters  would 
have  to  be  obtained  from  similarly  damaged  test  specimens.  Another 
possibility  would  be  to  represent  the  material  properties  as  a  function 
of  damage  level  and  incorporate  this  in  the  material  characterization. 
Of  course,  further  research  would  be  necessary  before  this  approach 
can  be  judged  feasible. 
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SECTION  6  STUDY  OF  LAMINATED  ELASTOPLASTIC  MATERIALS 


INTRODUCTION 

Thin  sheets  of  material  are  often  layered  (or  laminated)  together  to 
obtain  a  composite  which  exhibits  superior  performance.  While  elastic 
analyses  of  laminated  structures  have  been  thoroughly  studied,  only 
recently  have  elastic-plastic  materials  been  considered  [30,31  ].  The 
purpose  of  this  chapter  is  to  extend  the  previous  research  on  homogenous 
anisotropic  elastoplastic  materials  and  combine  it  with  the  incremental 
lamination  theory.  Then,  utilizing  the  previously  developed  computa¬ 
tional  tools,  the  complex  behavior  of  various  laminated  structures  will 
be  explored. 

In  this  section,  the  global  behavior  of  a  laminate  consisting  of 

isotropic  elastoplastic  layers  is  first  investigated  utilizing  the 
LAMINATE  program.  A  laminate  composed  of  orthotroic  layers  is  also 
analyzed  and  comparisons  drawn  to  the  isotropic  example.  The  ANPLAST 
finite  element  program  is  used  to  analyze  a  laminated  plate  with  a 
circular  hole,  for  both  isotropic  and  orthotropic  layer  construction. 

ELASTIC-PLASTIC  LAMINATE  BEHAVIOR:  ISOTROPIC  LAYERS 
The  global  response  of  a  laminate  consisting  of  isotropic  elastic- 
plastic  layers  was  studied.  Uniaxial  and  multiaxial  behavior  was  inves¬ 
tigated  using  LAMINATE. 
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Global  Response  Due  to  Uniaxial  Loading  and  Unloading 

The  uniaxial  response  of  a  laminate  consisting  of  two  isotropic 
elastic-plastic  materials  was  analyzed,  and  the  results  are  illustrated  in 
Fig.  6.1.  The  elastoplastic  material  properties  chosen  resemble  what 
one  might  find  for  a  typical  steel  and  aluminum  alloy.  The  work¬ 
hardening  behavior  is  assumed  to  be  of  the  Ramberg-Osgood  type.  Though 
the  analysis  requires  only  two  distinct  layers,  at  least  4  laminas 
would  be  required  for  symmetric  behavior. 

The  uniaxial  response  of  the  base  materials  was  also  included  to 
reinforce  the  concept  that  global  laminate  behavior  is  a  superposition 
of  individual  layer  response.  The  first  indication  of  this  averaging 
effect  is  the  two  distinct  break  points  for  the  laminate.  The  laminate 
first  experiences  yield  at  less  than  40,000  psi ■  This  behavior  is  due 
to  the  mismatch  in  elastic  materia?  properties.  Layer  #1  being  three 
times  stiffer  than  Layer  #2,  it  carries  a  significantly  greater  pro¬ 
portion  of  the  remote  laminate  loading.  Therefore,  it  approaches 
yield  at  a  much  greater  rate  than  does  a  specimen  consisting  of  a 
single  layer.  Unloading  is  also  included  to  demonstrate  that  permanent 
deformation  is  present  after  the  initial  yield.  It  is  interesting  to 
note  that  after  the  initial  yield,  the  slope  of  the  stress-strain 
response  is  nearly  constant.  This  is  due  to  the  dominance  of  elastic 
stiffness  contribution  of  Layer  #2. 

At  approximately  57000  psi.  Layer  #2  yields  and  the  entire 
laminate  becomes  elastic-plastic.  The  global  stress-strain  response 
now  becomes  non-linear  due  to  the  contribution  of  each  lamina's  non¬ 
linear  work-hardening  law.  It  is  apparent  that  the  laminate  response 


178 


NET  SECTION  STRESS  (PSD  (XIO 


1.000 


.800 


.600 


.HOO 


.200 


0.000 

0. 


Figure  6 


LAMINATE  #1 


0000 


.1000  .2000 
AXIAL  STRAIN 


.3000 


.HOOO 

(XIO  -1) 


.5000 


Layer  #1 

E  =  30  X  10^  psi 
=  50000  psi 

?  =  .  30 

n  =  9 


Laye'^ 

E  =  10  X  10^  psi 

"  =  50000  psi 

y 

:  ^  .  30 

n  =  7 


1  Uniaxial  load-deformation  response  of  an  isotropic 
laminate  (Laminate  #1). 
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tends  to  average  the  response  of  the  individual  layers.  The  laminate 
was  unloaded  again  to  demonstrate  the  resulting  permanent  deformation. 

The  effect  of  changing _the  layer  material  properties  on 
global  laminate  reponse  was  studied  next.  The  elastic  material 
properties  are  unchanged  from  the  previous  example,  only  the  uniaxial 
yield  stresses  for  each  material  was  changed.  Layer  #1  is  now  assumed 
to  have  a  yield  stress  twice  the  value  of  Layer  #2. 

This  change  in  material  properties  does  not  alter  the  general  behavior 
of  the  laminate,  as  illustrated  in  Fig.  6.2.  It  is  apparent  that  the 
elastoplastic  stress-strain  response  of  the  base  materials  are  more 
disperse,  yet  the  laminate  tends  to  average  the  responses.  The  initial 
yield  still  occurs  in  Layer  #1  (steel  alloy)  due  to  the  elastic 
disparity,  but  is  soon  followed  by  yielding  in  Layer  #2  (aluminum 
alloy).  Note  that  both  laminate  yield  points  occur  completely  between 
the  base  material  yields.  The  laminate  yield  points  for  the  previous 
example  stradled  the  individual  material  yields,  i.e.,one  higher,  one 
lower. 

The  ratio  of  axial  to  transverse  stresses  within  each  layer  is 
presented  in  Fig.  6.3.  The  presence  of  a  transverse  stress  in  each 
lamina  is  due  to  the  mismatch  of  the  contractural  strain  contribution 
in  each  layer.  During  the  elastic  loading  range,  the  T.ismatch  is 
directly  associated  to  the  difference  of  the  Poisson  ratios.  After 
Layer  #1  has  yielded,  this  mismatch  becomes  even  more  pronounced  due 
to  the  contribution  of  the  plastic  Poisson  ratio,  which  in  this  case 
is  ^  .  When  both  layers  have  become  elastic-plastic,  the  tendency 
is  for  the  transverse  stress  to  become  negligible  as  axial  loading  is 
continued. 
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Figure  6.2  Uniaxial  load-deformation  response  of  an  isotropic 
laminate  (Laminate  #2). 
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3  Axial  and  transverse  layer  stress  distribution  during 
uniaxial  loading  of  an  isotropic  laminate. 
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It  is  quite  clear  that  the  slope  of  the  plots  in  Fig.  6.3  are 
not  constant,  which  indicates  non-proportionality  of  the  normal 
stresses.  It  is  interesting  to  note  that  this  non-proportional  stress 
state  occurs  even  though  the  remote  laminate  loading  is  monotonic 
uniaxial.  This  non-proportionality  does  not  present  a  problem  since 
the  formulation  of  the  plasticity  flow  rule  is  entirely  incremental. 
These  results  suggest  that  a  deformation  theory  of  plasticity  would  not 
be  adequate  for  elastic-plastic  laminates. 

The  lamination  of  isotropic  elastic-plastic  layers  appears  to 
produce  a  composite  with  a  unique  response.  The  global  behavior  of 
the  laminate  is  a  function  of  the  elastic  and  post  yield  plasticity 
parameters  of  the  individual  laminas.  It  has  been  demonstrated  that 
many  different  laminate  behaviors  can  be  characterized  with  knowledge 
of  the  base  layer  material  properties  and  use  of  the  LAMINATE  program. 
Ideally,  one  could  combine  several  materials  and  arrive  at  a  composite 
behavior  which  envelopes  the  best  attributes  of  each  base  material. 

Laminate  Response  Due  to  Cyclic  Loading 

The  global  response  of  a  laminate  under  a  uniaxial  load  cycled 
between  fixed  load  limits  was  investigated.  Laminate  #1  elastoplastic 
properties  were  selected  for  the  analysis.  Fig.  6.4  illustrates  the 
uniaxial  laminate  response  when  subjected  to  cyclic  loading.  After 
both  layers  have  yielded  (A,  B)  and  considerable  work-hardening  has 
occurred,  the  laminate  was  unloaded  (C)  and  cycled  three  times  between 
fixed  limits  which  are  equal  in  tension  and  compression. 
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Figure  6.4  Stress-strain  response  of  an  isotropic  laminate 
cycled  between  fixed  load  limits. 
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The  unloading  process  from  point  C  to  D  proceeds  in  an  elastic 
manner  as  is  consistent  with  elastoplastic  theory  but  yields  at  point 
D,  rather  than  continuing  elastically  until  the  lower  load  limit  is 
reached.  This  appears  to  be  a  contradiction  of  the  isotropic  hardening 
assumption  which  predicts  self-similar  expansion  of  the  yield  surface. 
This  behavior  can  be  accounted  for  if  one  recalls  that  each 
lamina  is  constrained  by  its  own  individual  elastic-plastic  yield 
criterion  and  isotropic  work-hardening  relationships.  The  interelation 
of  these  different  elastoplastic  materials  through  the  lamination 
assumptions  results  in  this  complex  behavior. 

During  the  second  and  third  cycles,  the  hysteresis  loop  appears  to 
become  narrower,  and  the  corresponding  tension  and  compression  yield 
points  move  closer  to  the  fixed  limits.  If  the  cyclic 
loading  continues,  it  has  been  demonstrated  that  the  hysteresis  loop 
eventually  collapses  to  a  straight  line  bounded  by  the  fixed  load 
limits. 

The  degree  of  permanent  deformation  exhibited  by  the  individual 
layers  cannot  be  discerned.  The  layer  stresses  are  plotted  against 
net  section  strain  in  Fig.  6.5.  It  is  clear  that  Layer  #1  underwent 
elastic-plastic  work-hardening,  while  Layer  #2,  the  aluminum  alloy, 
stayed  entirely  elastic.  The  explanation  of  this  phenomenon  lies  in 
the  much  greater  elastic  stiffness  of  Layer  #1.  As  a  result  of  the 
first  order  laminate  assumptions,  stiff  layers  tend  to  stress  to  a 
greater  extent  than  more  flexible  layers.  An  interesting  aspect  of 
this  behavior  is  that  the  elastic-plastic  response  is  dominated  by 
the  relative  magnitude  of  the  elastic  properties,  and  not  just  the 
plasticity  parameters. 
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LAMINATE  #1 
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Layer  stress-strain  response  of  an  isotropic  laminate 
cycled  between  fixed  load  limits. 


A  Study  of  the  Laminate  Yield  Surface 

The  results  from  the  previous  section  infer  that  Laminate  #1  does 
not  exhibit  isotropic  hardening.  Rather,  a  form  of  kinematic 
hardening  [ 8 ]  or  a  combination  of  isotropic  and  kinematic  hardening 
behavior  seems  apparent.  Figure  6,6  depicts  the  translation  and 
expansion  of  the  laminate  yield  surface  during  the  cyclic  loading 
represented  in  Fig.  6.4,  The  yield  surface  is  defined  uy  utilizing 
the  multiaxial  capability  of  the  LAMINATE  program.  The  first  plot 
(  O)  is  for  the  initial  yield  surface  iA),  which  '>s  centered  about 
the  origin.  The  next  plot  (A)  is  a  representation  of  the  laminate 
yield  surface  immediately  after  unloading  (C).  The  translation  of  the 
yield  locus  accounts  for  the  different  values  of  yield  in  tension  and 
compression  (C,  D).  It  is  also  clear  that  the  yield  surface  is 
expanding,  signifying  a  combined  isotropic-kinematic  work-hardening. 

The  third  yield  locus  (  +  )  is  representative  of  the  laminate  after 
re-load  (E).  The  yield  surface  has  translated  towards  negative 

with  ve^'y  little  expansion.  The  last  plot  depicts  the  laminate 
yield  surface  after  many  cycles,  when  the  hysteresis  loop  has  collapsed 
to  a  line.  It  is  apparent  that  the  anisotropy  of  yield  induced  by  the 
initial  load  cycles  eventually  dissipates  as  the  laminate  continues  to 
work-harden. 

A  multiaxial  loading  path  and  subsequent  work-;.  euing  was 
considered  next.  A  set  of  yield  surfaces  at  three  consecutive  levels 
of  load  is  illustrated  in  Fig.  6.7.  The  initial  yield  surface  takes 
the  familiar  form  of  an  ellipsoid  centered  about  the  origin  (O).  The 
yield  surfaces  at  points  (A)  and  (+)  indicate  both  kinematic  and 
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Figure  6.6  Illustration  of  the  yield  surface  for  an  isotropic 
laminate  during  uniaxial  cyclic  loading  between 
fixed  load  limits. 


188 


aimnate  Load  Path 


Figure  6.7  Illustration  of  the  yield  surface  for  an  isotropic 
laminate  during  a  multiaxial  non-proportional 
loading  path. 


isotropic  hardening,  which  is  consistent  with  the  cyclic  loading 
results. 

Though  only  isotropic  hardening  behavior  is  assumed  for  each 
individual  layer,  the  laminate  demonstrates  a  more  coniplex  response. 

In  fact,  depending  on  the  particular  make-up  of  a  laminate,  it  is 
likely  that  global  load- deformation  response  cannot  be  described  by 
classical  elastoplasticity  theory.  Fortunately,  the  researcher  need 
not  bother  with  the  laborious  task  of  characterizing  the  laminate  as  a 
homogenous  material;  rather,  one  need  only  have  a  knowledge  of  the 
simpler  lamina  properties,  and  a  properly  formulated  lamination  theory. 

Path  Dependency  Due  to  Multiaxial  Loading 

Two  sets  of  complex  loading  paths  involving  multiaxial  stresses 
were  studied.  Figure  6.8  illustrates  the  normal  and  shear  strain 
dependence  on  load  path  selection.  These  plots  draw  a  strong 
similarity  to  the  study  conducted  in  Section  ^  which  concerned  strain 
dependence  on  load  path  for  a  homogenous  elastic-plastic  material 
(Fig.  2.13).  The  same  conclusions  can  be  reached  for  the  laminate  as 
for  the  homogenous  material;  stress-strain  response  is  not  unique  for 
complex  loading  paths. 

ELASTIC-PLASTIC  LAMINATE  BEHAVIOR:  ORTHOTROPIC  LAYERS 
Laminates  consisting  of  orthotropic  elastoplastic  layers  were 
considered  and  comparisons  made  with  isotropic  laminate  behavior.  This 
topic  is  of  practical  interest  since  many  metal  matrix  composite 
applications  require  lamination. 
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Figure  6.8  Deformation  response  of  an  isotropic  laminate 
due  to  different  complex  loading  paths. 
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Global  Response  Due  to  Uniaxial  Loading  and  Unloading 

The  uniaxial  response  of  a  laminate  [0°/90°]g  consisting  of  two 
orthotropic  elastic-plastic  materials  (B-Al)  is  pictured  in  Fig.  6.9. 
The  laminate  exhibits  an  initial  linear  behavior,  then  becomes  non¬ 
linear  after  the  90°  ply  yields.  Though  the  0°  ply  eventually  yields 
its  stiffness  contribution  remains  essentially  elastic  due  to  the 
fiber-load  alignment.  Note  how  the  laminate  response  again  falls 
between  the  homogenous  layer  material  behavior. 

The  unloading  behavior  initially  demonstrates  the  classical 
elastic  response,  but  at  approximately  20,000  psi  re-yields.  It  is 
clear  that  the  yield  surface  has  been  distorted  during  the  tensile 
work  hardening,  introducing  some  degree  of  additional  anisotropy. 

After  the  load  is  completely  removed,  a  permanent  net  section  laminate 
strain  is  incurred,  substantiating  the  idea  that  plastic  deformation 
is  occurring  in  the  various  layers.  This  type  of  behavior  has  been 
observed  by  Dvorak  [30]  and  Sova  [32]. 

The  stresses  in  the  laminas  during  laminate  deformation  are 
illustrated  in  Fig.  6.10.  Transverse  stresses  in  the  individual 
layers  are  apparent  as  was  the  case  with  the  isotropic  laminate.  The 
contractural  mismatch  is  due  to  the  orthotropic  nature  of  each  layer, 
and  the  particular  lay-up  design  [0°/90°]j.  The  ratio  of  the  principal 
stress  components  denoted  by  the  slopes  of  each  plot  is  obviously  not 
constant.  This  again  reinforces  the  non-proportionality  of  the  stress 
state. 

The  induction  of  transverse  stress  into  the  laminas  during 
laminate  deformation  has  an  effect  on  subsequent  yield.  For  example. 
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Figure  6.9  Uniaxial  load-deformation  response  of  an 
orthotropic  B-Al  laminate  [0°/90°]g. 


93 


Figure  6.10  Axial  and  transverse  layer  stress  distribution 
during  monotonic  uniaxial  loading  of  the  B*A1 
laminate. 


the  0°  ply  has  the  larger  axial  stress,  but,  due  to  the  anisotropy  of 
the  yield  criterion,  this  stress  component  does  not  significantly 
contribute  to  yielding.  Conversely,  the  transverse  stress  present  in 
this  ply  does  contribute. 

Laminate  Response  Due  to  Cyclic  Loading 

The  global  response  of  this  same  laminate  subjected  to  a  uniaxial 
load  cycled  between  fixed  load  limits  was  studied.  The  B-Al  results 
are  illustrated  in  Fig.  6.11.  After  initial  yield,  the  laminate 
work-hardened  until  unloading  at  point  B.  The  load  was  removed  and  the 
laminate  underwent  an  equal  amount  of  compression.  The  laminate  was 
then  reloaded  until  the  upper  limit  was  reached.  The  process  was 
repeated  for  another  two  cycles. 

The  results  are  very  similar  to  those  for  the  isotropic  laminate. 
Unloading  initially  occurs  elastically,  but  reverse  yielding  occurs 
before  the  bottom  load  limit  is  reached,  again  indicating  that  non¬ 
isotropic  hardening  occurs.  The  hysteresis  loops  for  the  Boron- 
Aluminum  laminate  are  initially  narrower  and  collapse  to  a  line  faster 
than  for  the  isotropic  laminate.  This  behavior  is  associated  with  the 
greater  ductility  demonstrated  by  the  isotropic  materials.  Due  to  the 
fiber  orientation  of  the  0°  ply, plastic  deformation  is  retarded. 

The  stress-strain  behavior  in  the  individual  laminas  are  plotted 
in  Fig.  6.12.  The  90°  ply,  which  yields  first,  experiences  considerable 
plastic  deformation  as  suggested  by  the  large  hysteresis  loops.  The 

0°  ply  demonstrates  almost  purely  elastic  behavior  throughout.  This 
behavior  is  consistent  with  the  view  that  very  little  permanent 
deformation  can  occur  parallel  to  the  stiff  fibers. 
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Figure  6.11  Stress-strain  response  of  the  B-Al  laminate 
cycled  between  fixed  load  limits. 
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90°  ply 


Figure  6.12  Layer  stress-strain  response  of  the  B-Al  laminate 
cycled  between  fixed  load  limits. 
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It  is  interesting  to  note  that  while  the  isotropic  laminate 
behavior  was  dominated  by  the  ratio  of  the  elastic  stiffnesses,  the 
Boron-Aluminum  laminate  is  not  strongly  dependent  on  elastic  properties. 
Instead,  the  Boron-Aluminum  is  dominated  by  the  anisotropy  of  yield, 
and  the  corresponding  orientation  of  the  individual  laminas. 

A  Study  of  the  Laminate  Yield  Surface 

The  laminate  unloading  and  cyclic  loading  behavior  demonstrates 
that  the  work-hardening  is  not  isotropic.  Figure  6.13  depicts  four 
yield  surfaces  at  various  load  points  from  the  cyclic  loading  example 
in  Fig.  6.11. 

The  initial  yield  locus  (©)  in  Fig.  6.13  is  almost  square  in 
shape.  This  is  quite  different  than  that  depicted  for  the  isotropic 
laminate.  This  contrasting  behavior  is  due  to  the  fibers  being 
aligned  with  each  principal  axis.  This  results  in  biaxial  stresses 
having  a  negligible  effect  on  the  laminate  yield  surface.  Only  the 
stress  within  each  layer  that  is  transverse  to  the  fibers  contributes 
to  the  yield.  The  yield  surface  reported  by  Dvorak  [30]  takes  a 
similar  form. 

The  second  plot  {A)  represents  the  laminate  yield  surface  after 
tensile  work-hardening  to  load  point  B.  The  yield  surface  has  both 
translated  in  an  approximate  rigid  body  sense,  and  experienced 
distortion.  The  third  locus  (+)  represents  the  yield  surface  at  the 
bottom  of  the  first  load  cycle  (D)  and  has  also  undergone  a  translation 
and  further  distortion.  Finally,  at  the  end  of  the  third  loading 
cycle,  the  yield  locus  (x)  extends  from  the  maximum  to  the  minimum 
load  limit. 
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Figure  6.13  Illustration  of  the  yield  surface  for  the  B-Al  laminate  during 
uniaxial  cyclic  loading  between  fixed  load  limits. 


It  is  obvious  that  B-Al  laminate  work-hardening  is  highly 
anisotropic.  The  anisotropy  of  the  work-hardening  is  not  only  induced 
by  the  rigid  body  translation,  but  the  distortion  of  the  yield  locus 
as  well.  Note  that  the  yield  surface  tends  to  distort  in  the  direction 
of  the  applied  load,  but  expands  a  relatively  small  amount  in  the 
transverse  direction.  It  is  likely  that  this  type  of  deformation  is 
due  to  the  relatively  stiff  behavior  of  the  fibers  in  resisting 
transverse  plastic  deformation. 

The  yield  and  work-hardening  behavior  studied  in  this  section 
provides  justification  for  not  analyzing  the  laminate  on  the 
macroscopic  scale.  The  theoretical  mechanics  necessary  to  describe 
this  unusual  yield  surface  and  subsequent  work-hardening  would  be  much 
too  cumbersome.  Understanding  the  elastic-plastic  behavior  at  the 
“mini -level"  and  superimposing  the  responses  via  a  lamination  theory, 
appears  to  offer  the  straightforward  approach. 

Path  Dependency  Due  to  Multiaxial  Loading 

Load  path  dependency  of  strain  was  investigated  to  demonstrate 
the  non-unique  stress-stain  response.  Two  complex  loading  paths 
involving  sequenced  multiaxial  loads  were  applied  to  the  B-Al  laminate. 
Figure  6.14  illustrates  the  normal  and  shear  strain  dependence  on  the 
particular  load  path.  These  load  paths  are  similar  to  those  for  the 
isotropic  laminate,  in  that  there  is  a  wide  disparity  in  the  final 
deformation  state. 

The  similarity  of  the  isotropic  homogenous,  isotropic  laminate, 
and  anisotropic  B-Al  laminate  response  confirms  that  the  non¬ 
proportionality  of  the  stresses  has  a  dominant  effect  on  the  final 
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Figure  6.14  Deformation  response  of  the  B-Al  laminate  due  to 
different  complex  loading  paths. 
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deformation  state,  regardless  of  anisotropy  or  material  construction. 


ANALYSIS  OF  A  LAMINATED  SHEET  WITH  A  CIRCULAR  HOLE:  ISOTROPIC  LAYERS 
The  ANPLAST  finite  element  program  was  used  to  analyze  a  laminated 
sheet  with  a  circular  hole.  The  laminated  shee  is  constructed  of  two 
isotropic  layers,  as  characterized  in  Fig.  6.1.  The  geometry  and 
loading  of  the  finite  specimen  is  illustrated  in  Fig.  4.1. 

A  Study  of  Plastic  Zone  Growth  Within  Each  Layer 

Growth  of  the  plastic  zone  within  each  layer  is  illustrated  in 
Fig.  6.15.  As  expected,  yielding  is  initiated  in  Layer  #1  (Steel 
Alloy)  due  to  its  much  greater  elastic  stiffness.  The  initial  yield 
is  located  at  the  perimeter  of  the  hole,  along  the  symmetry  line, 
where  the  hoop  stress  (oyy)  a  maximum.  Yield  is  initiated  in  Layer 
#2  (Aluminum  Alloy)  at  a  greater  remote  load. 

The  shape  and  growth  of  the  plastic  zone  within  each  layer 
parallels  the  behavior  of  a  homogenous  isotropic  material,  as  was 
analyzed  in  Section  4.  It  is  apparent  that  the  two  layers  interact 
in  such  a  way  that  the  character  of  the  plasticity  is  a  function  of 
the  individual  layer  properties,  but  the  magnitude  of  plastic  zone 
growth  is  a  superposition  of  the  two  responses.  On  the  element  level 
total  strain  within  each  layer  is  the  same,  yet  different  levels  of 
plastic  zone  growth  are  found.  This  indicates  that  the  elastic  (e®) 
and  plastic  (e^)  components  of  strain  can  be  different  for  each  layer, 
but  their  sum  must  be  the  same. 

A  more  quantitative  presentation  of  plastic  zone  growth  is 
obtained  by  plotting  contours  of  effective  stress  as  shown  in  Fig. 
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Figure  6.15  Plastic  zone  growth  in  each  layer  of  an 
isotropic  laminated  sheet  subjected  to 
elastic-plastic  loading. 


203 


Layer  #2 

o.  =  33500  = 


Layer  #1 


Figure  6.15  continued 
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6.16.  These  contours  reiterate  the  different  states  of  plastic 
deformation  in  each  layer. 

A  Study  of  Residual  Stresses  and  Strains  Due  to  Unloading 

The  specimen  was  first  unloaded  by  taking  a  small  negative  loading 
step,  which  allowed  the  elements  to  become  elastic.  The  small  unloading 
increment  also  provided  an  opportunity  for  re-assembly  of  the  stiffness 
matrix.  The  load  was  removed  in  ten  equal  increments.  In  this 
particular  case,  both  layers  and  all  elements  became  elastic  during  the 
unloading,  and  did  not  re-yield. 

The  residual  stress  (o^y)  distributions  are  illustrated  in  Fig. 

6.17.  Layer#l  exhibits  the  expected  compressive  residual  stress  component, 
while  the  magnitude  of  stress  in  Layer  #2  is  considerably  less.  The 
reason  for  this  behavior  involves  the  greater  degree  of  permanent 
deformation  experienced  by  Layer  #1.  In  addition,  since  Layer  #1  is 
stiffer,  it  will  unload  faster  for  an  equal  increment  of  laminate 
strain.  It  is  interesting  to  hypothesize  that  while  the  compressive 
residual  stress  in  Layer  #1  may  be  sufficient  to  retard  crack 
initiation  due  to  an  overload  this  may  not  be  the  case  for  Layer  #2. 

ANALYSIS  OF  A  LAMINATED  SHEET  WITH  A  CIRCULAR  HOLE:  ORTHOTROPIC  LAYERS 
An  analysis  of  the  laminated  sheet  with  a  circular  hole 
subjected  to  uniaxial  remote  loading  was  conducted  with  ANPLAST.  The 
laminated  sheet  consisted  of  two  layers  of  B-Al  [0°/90°]^. 

A  Study  of  Plastic  Zone  Growth  Within  Each  Layer 

Growth  of  the  plastic  zone  within  each  layer  is  illustrated  in 
Fig.  6.18.  The  initiation  of  permanent  deformation  occurs  in  the  90° 
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Effective  stress  contours  in  each  layer  for  the  isotropic  laminated 
sheet  subjected  to  elastic-plastic  loading. 
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Figure  6.17  ANPLAST  residual  stresses  in  each  layer  for  the 
isotropic  laminated  sheet  with  a  circular  hole 
subjected  to  elastic  unloading. 
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Figure  6.18  Plastic  zone  growth  in  each  ply  for  the  B-Al 
laminated  sheet  [0°/90°]^  with  a  circular 
hole  subjected  to  elastic*plastic  loading. 
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Figure 


.18  continued 
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layer.  Unlike  the  isotropic  laminate  results,  the  initiation  of 
yielding  is  not  caused  by  the  ratio  of  the  layer  elastic  stiffness. 
Rather,  the  very  strong  plastic  anisotropy  and  ply  orientation  dominate. 

The  shape  and  growth  of  the  plastic  zones  within  each  layer 
parallels  the  response  observed  in  the  orthotropic  homogenous  problem. 
The  plastic  zone  in  the  90°  ply  tends  to  grow  across  the  width  of  the 
specimen,  while  the  zone  in  the  0°  ply  extends  in  a  more  vertical 
pattern.  Note  that  net  section  yielding  is  prevalent  in  the  90°  ply  at 
a  remote  load  of  3^000  psi,  while  at  the  same  time  the  plastic  zone  in 
the  0°  ply  remains  constrained  to  the  vicinity  of  the  hole.  Though 
the  plastic  zones  in  each  layer  take  their  characteristic  shape 
observed  in  Section  4,  the  rate  of  progression  within  each  layer  is 
controlled  by  the  superposition  of  the  two  responses. 

The  contour  of  effective  stress  in  Fig.  6.19  clearly  demonstrates 
the  much  greater  magnitude  of  plastic  deformation  present  in  the  90° 
ply.  For  the  elastic  loading  increment,  the  contour  of  the  0°  ply 
indicates  that  the  initiation  of  yield  did  not  occur  at  the  symmetry 
line,  but  along  the  hole  boundary  at  approximately  e  =  10°.  This  is 
consistent  with  the  behavior  observed  in  the  analysis  of  the 
homogenous  orthotropic  material. 

Dvorak  [33]  has  also  analyzed  a  laminated  plate  [0°/90°]^  with  a 
circular  hole  for  a  FP-aluminum  metal  matrix  composite.  The  plastic 
zone  visualization  for  this  analysis  is  illustrated  in  Fig.  6.20.  Note 
the  similarity  of  plastic  zone  shape,  especially  at  the  lower  load. 

It  is  interesting  to  note  that  Dvorak  models  the  composite  as  a 
micromechanical  structure,  in  contrast  to  the  simpler  macroscopic 
approach  taken  in  this  research. 


210 


ELASTIC  ‘  l^OOO  MAX.  LOAD  RESIDUAL 


6.19  Effective  stress  contours  in  each  ply  for  the  B-Al  laminated 
sheet  subjected  to  elastic-plastic  loading. 
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Figure  6.20  Plastic  zone  growth  calculated  by  Dvorak  [33]  for 
an  orthotropic  FP-Al  laminated  sheet  [0°/90‘']  with 
a  circular  hole  subjected  to  elastic-plastic 
loading. 
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A  Study  of  Residual  Stresses  and  Strains  Due  to  Unloading 

The  B-Al  laminate  was  unloaded  in  the  same  fashion  as  the  isotropic 
laminate,  a  small  unloading  step  followed  by  ten  equal  unloading 
increments.  After  the  first  unloading  step,  both  layers  became 
elastic  and  remained  elastic  until  the  end  of  the  last  loading 
increment  when  4  elements,  in  proximity  to  the  hole  boundary, 
re-yielded. 

The  residual  stress  distributions  are  illustrated  in  Fig. 

6.21.  The  typical  compressive  residual  stress  exists  in  the  90°  ply, 
but  a  residual  tensile  stress  is  present  in  the  0°  ply.  Since  is 
also  tensile  in  the  0°  ply,  these  results  suggest  that  crack  initiation 
due  to  cyclic  loading  may  be  prevalent  in  this  ply. 

Experience  and  intuition  are  not  sufficient  in  analyzing  these 
complex  elastic-plastic  laminates.  As  the  geometries  and  loadings 
become  even  more  complex,  the  need  for  a  practical  computational  tool, 
such  as  ANPLAST,  becomes  obvious. 
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Figure  6.21  ANPLAST  residual  stresses  in  each  ply  for 
the  B-Al  laminated  sheet  with  a  circular 
hole  subjected  to  elastic  unloading. 
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SECTION  7 


CONCLUSIONS 


The  purpose  of  this  research  was  to  develop  a  numerical  tool  to 
investigate  the  role  of  anisotropy  in  the  elastic-plastic  behavior 
of  structural  components. 

The  plane  stress  incremental  orthotropic  plasticity  relations  for 
macroscopic  deformation  of  thin  sheets  have  been  properly  formulated.  The 
plasticity  coefficients  *^66^’  effective  yield  stress  (o^) 

and  particular  work-hardening  law  (e*^  vs  o)  control  the  shape  and 
growth  of  the  yield  surface  as  well  as  the  associated  plastic  flow. 

The  incremental  formulation  provides  the  necessary  flexibility  to 
analyze  nonproportional  loading,  including  unloading. 

The  orthotropic  plasticity  flow  rule  was  formulated  so  that  it 
is  compatible  with  the  finite  element  method  and  directly  implemented 
into  the  program  ANPLAST.  Both  the  constant  strain  triangle  and 
isoparametric  quadralateral  elements  provide  satisfactory  results. 

The  non-linear  problem  is  solved  in  an  incremental  manner  which  allows 
complete  non-proportionality  of  loading.  The  pre-processing  and 
post-processing  software  provides  an  efficient  compliment  to  the 
overall  ANPLAST  analysis  capabilities. 

The  analysis  of  the  sheet  with  a  circular  hole  provides  a  bench¬ 
mark  for  the  elastic  and  isotropic  elastic-plastic  capabilities  of 
ANPLAST.  The  isotropic  and  orthotropic  elastic  capability  of  ANPLAST 
demonstrates  excellent  agreement  with  the  theoretical  solution.  The 
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isotropic  elastoplastic  performance  of  ANPLAST  compares  well  with 
available  theoretical  and  experimental  solutions.  Close  concordance 
with  the  experimental  results  of  Rizzi  [24]  for  a  strongly  anisotropic 
B-Al  metal -matrix  composite  verifies  the  adequacy  of  the  anisotropic 
formulation  and  demonstrates  the  practicality  of  applying  ANPLAST. 

The  elastic  analysis  of  the  cracked  panel  illustrates  the  conven¬ 
ience  of  utilizing  an  energy  method  for  stress  intensity  calibration 
of  an  anisotropic  material.  An  elastoplastic  analysis  reveals  that 
plastic  flow  in  the  vicinity  of  the  crack  is  strongly  influenced  by 
the  amount  of  anisotropy  and  the  orientation  of  the  principal  material 
direction.  This  suggests  that  the  mechanisms  which  control  crack 
propagation  are  firmly  linked  to  material  directionality  and  the 
resulting  plasticity. 

A  number  of  strong  similarities  exist  between  the  elastoplastic 
results  for  the  hole  and  crack.  The  progression  of  the  plastic  zone 
and  related  plastic  flow  is  controlled  by  the  same  parameters.  In 
both  cases  a  strong  compressive  longitudinal  residual  stress  is  evident 
at  the  perimeter  of  each  discontinuity,  regardless  of  the  degree  of 
anisotropy. 

The  response  of  elastoplastic  laminated  materials  exhibits 
extremely  complex  behavior.  The  study  of  both  isotropic  and  orthotro¬ 
pic  layers  leads  to  the  conclusion  that  global  laminate  response  is  a  *• 

non-linear  superposition  of  individual  elastic  and  elastic-plastic 
lamina  properties.  A  laminate  consisting  of  orthotropic  layers 
exhibits  an  initial  yield  surface  which  is  quite  different  from  that 
of  each  lamina.  Loading  of  either  laminate  construction  beyond  its 
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initial  yield  results  in  a  translation  and  distortion  of  the  yield 
surface.  When  the  laminate  is  loaded  between  fixed  load  limits,  the 
yield  surface  eventually  stabilizes  after  a  number  of  load  cycles. 

It  is  clear  that  it  is  difficult  to  apply  classical  elastoplasticity 
theory  to  the  global  deformation  of  laminates.  But  it  is  shown  that 
the  most  straightforward  approach  to  the  solution  of  these  problems  is 
the  characterization  of  elastic-plastic  behavior  at  the  mini-level 
and  superposition  of  the  lamina  responses  via  a  legHimate  lamination 
theory. 

The  analysis  of  the  laminated  sheet  with  a  circular  hole  exempli¬ 
fies  the  usefulness  of  the  implementation  of  the  incremental  lamination 
theory  into  ANPLAST.  The  results  demonstrate  the  significantly 
different  plastic  deformation  occurring  within  each  layer. 

Anisotropy  can  have  a  very  complex  influence  on  the  elastoplastic 
response  of  structures.  But  with  aid  of  a  numerical  tool  such  as 
ANPLAST  they  can  be  made  tractable  and  their  interesting  aspects 
conveniently  elucidated. 
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